1 DATA

## Longitudinal data
data_sum <- loadlongitudinaldata(dataset = "DATA_Adults_G1G29.csv", rm_generation1 = 1,rm_generation2 = 7,rm_generation3 = 29)

## Phenotyping steps
data_G0 <- loadfitnessdata(dataset = "Selection_Phenotypage_G0_G7_G8.csv", generation = "G1")
data_G7 <- loadfitnessdata(dataset = "Selection_Phenotypage_G0_G7_G8.csv", generation = "G7")
data_G29 <- loadfitnessdata(dataset = 
                            "PERFORMANCE_Comptage_adultes_G13G14G15G16G17G18G19G20G21G22G23G24G25G26G27G28G29.csv",
                            generation = "29")
head(data_sum)
##   Line Fruit_s Generation          Step  N Nb_adults       sd    fitness
## 1  CE1  Cherry          2 first_prepool 20   9.15000 6.123939 -0.7819784
## 2  CE1  Cherry          3 first_prepool 10  15.30000 7.631077 -0.2678794
## 3  CE1  Cherry          4 first_prepool  8  15.00000 7.782765 -0.2876821
## 4  CE1  Cherry          5 first_prepool  6  14.50000 6.284903 -0.3215836
## 5  CE2  Cherry          2 first_prepool 20   7.75000 7.758696 -0.9480394
## 6  CE2  Cherry          3 first_prepool  7  14.42857 8.303757 -0.3265219
##   se_fitness
## 1  0.1496562
## 2  0.1577228
## 3  0.1834415
## 4  0.1769518
## 5  0.2238577
## 6  0.2175215
head(data_G0)
##     Treatment Line Fruit_s Nb_eggs Nb_adults SA Emergence_rate
## 993    Cherry  CE1      GF      76         6  0     0.07894737
## 994    Cherry  CE1      GF      89        17  0     0.19101124
## 995    Cherry  CE1      GF      57        12  0     0.21052632
## 996    Cherry  CE1      GF     172        24  0     0.13953488
## 997    Cherry  CE1      GF     173        33  0     0.19075145
## 998    Cherry  CE1      GF      91        18  0     0.19780220
head(data_G7)
##    Treatment Line    Fruit_s Nb_eggs Nb_adults SA Emergence_rate
## 3 Strawberry  CR4  Cranberry     152        68  0      0.4473684
## 4  Cranberry  CR4  Cranberry     246        25  1      0.1016260
## 5     Cherry  CR4  Cranberry     238        29  0      0.1218487
## 6     Cherry  CR4  Cranberry     166        23  0      0.1385542
## 8  Cranberry  FR3 Strawberry     204         5  0      0.0245098
## 9 Strawberry  FR3 Strawberry     124        45  1      0.3629032
head(data_G29)
##       Treatment Line Fruit_s Nb_eggs Nb_adults SA Emergence_rate
## 5392 Strawberry  CEA  Cherry     196        16  0     0.08163265
## 5393 Strawberry  CEA  Cherry     192        30  0     0.15625000
## 5394 Strawberry  CEA  Cherry     160        17  0     0.10625000
## 5395 Strawberry  CEA  Cherry     106         9  0     0.08490566
## 5396 Strawberry  CEA  Cherry     119        14  0     0.11764706
## 5397 Strawberry  CEA  Cherry     204        24  0     0.11764706
dim(data_G29)
## [1] 990   7
#Compute log change
data_logchange <- computelogchange(fitness_dataset_intermediate = data_G7, fitness_dataset_final = data_G29)
#Lines: CE2 and CR2 do not pass Geary's test (the threshold of a standardized mean greater than 3)
data_logchange$Line_type  <- ifelse(data_logchange$Line == "CR2"| data_logchange$Line == "CE2", "solid","dashed")


#Formatting for testing correlation
TEMP_dataG7_CheCran <- formattinglogchange(data_logchange, "7", "Cherry", "Cranberry")
TEMP_dataG7_CranStraw <- formattinglogchange(data_logchange, "7", "Cranberry", "Strawberry")
TEMP_dataG7_StrawChe <- formattinglogchange(data_logchange, "7", "Strawberry", "Cherry")

TEMP_dataG29_CheCran <- formattinglogchange(data_logchange, "29", "Cherry", "Cranberry")
TEMP_dataG29_CranStraw <- formattinglogchange(data_logchange, "29", "Cranberry", "Strawberry")
TEMP_dataG29_StrawChe <- formattinglogchange(data_logchange, "29", "Strawberry", "Cherry")

2 ADAPTATION

2.1 Analysis

######## Models
mod_nul <- lme4::lmer(fitness ~ 1 + (1|Generation), 
            weights = N, data = data_sum)

mod_Fruits <- lme4::lmer(fitness ~ Fruit_s + (1|Generation), 
            weights = N, data = data_sum)

mod_Generation <- lme4::lmer(fitness ~ Generation + (1|Generation), 
            weights = N, data = data_sum)

mod_Step <- lme4::lmer(fitness~ Step + (1|Generation), 
            weights = N, data = data_sum)

mod_Lines <- lme4::lmer(fitness~ Line + (1|Generation), 
            weights = N, data = data_sum)

mod_Generation_Step <- lme4::lmer(fitness~ Step + Generation + (1|Generation),  
            weights = N, data = data_sum)

mod_Generation_Fruits <- lme4::lmer(fitness ~ Generation + Fruit_s + (1|Generation), 
            weights = N, data = data_sum)

mod_Fruits_Step <- lme4::lmer(fitness ~ Step + Fruit_s + (1|Generation), 
            weights = N, data = data_sum)


mod_Generation_Fruits_Step <- lme4::lmer(fitness ~ Generation + Step + Fruit_s + (1|Generation), 
                              weights = N, data = data_sum)


mod_Doubleinteraction <- lme4::lmer(fitness ~ Generation + Step + Fruit_s + 
                              Fruit_s:Generation + 
                              (1|Generation),  
                              weights = N, data = data_sum)

mod_Double_gene <- lme4::lmer(fitness ~ Generation*Fruit_s + 
                              (1|Generation), 
                              weights = N, data = data_sum)

mod_Double <- lme4::lmer(fitness ~ Step*Fruit_s + 
                              (1|Generation), 
                              weights = N, data = data_sum)

mod_Double_Generation <- lme4::lmer(fitness ~ Generation + Step*Fruit_s + 
                              (1|Generation), 
                              weights = N, data = data_sum)

mod_Tripleinteraction <- lme4::lmer(fitness ~ Generation*Step*Fruit_s + 
                              (1|Generation), 
                              weights = N, data = data_sum)

mod_DoublePopulationGeneration <- lme4::lmer(fitness ~ Generation*Line + 
                              (1|Generation), 
                              weights = N, data = data_sum)



# AICC and Loglik
MuMIn::AICc(mod_nul,mod_Fruits,mod_Generation,mod_Step,mod_Lines,mod_Generation_Step,mod_Generation_Fruits,mod_Fruits_Step,mod_Generation_Fruits_Step,mod_Double_gene,mod_Doubleinteraction,mod_Double,mod_Double_Generation,mod_Tripleinteraction,mod_DoublePopulationGeneration)
##                                df     AICc
## mod_nul                         3 339.7756
## mod_Fruits                      5 333.0398
## mod_Generation                  4 334.4473
## mod_Step                        5 314.8770
## mod_Lines                      31 366.1793
## mod_Generation_Step             6 322.8487
## mod_Generation_Fruits           6 328.4534
## mod_Fruits_Step                 7 309.1298
## mod_Generation_Fruits_Step      8 317.0215
## mod_Double_gene                 8 342.2655
## mod_Doubleinteraction          10 330.9505
## mod_Double                     11 317.2153
## mod_Double_Generation          12 325.1311
## mod_Tripleinteraction          20 354.5837
## mod_DoublePopulationGeneration 60 504.0109
anova(mod_nul,mod_Fruits,mod_Generation,mod_Step,mod_Lines,mod_Generation_Step,mod_Generation_Fruits,mod_Fruits_Step,mod_Generation_Fruits_Step,mod_Doubleinteraction,mod_Double,mod_Double_Generation,mod_Tripleinteraction,mod_DoublePopulationGeneration,mod_Double_gene)
## Data: data_sum
## Models:
## mod_nul: fitness ~ 1 + (1 | Generation)
## mod_Generation: fitness ~ Generation + (1 | Generation)
## mod_Fruits: fitness ~ Fruit_s + (1 | Generation)
## mod_Step: fitness ~ Step + (1 | Generation)
## mod_Generation_Step: fitness ~ Step + Generation + (1 | Generation)
## mod_Generation_Fruits: fitness ~ Generation + Fruit_s + (1 | Generation)
## mod_Fruits_Step: fitness ~ Step + Fruit_s + (1 | Generation)
## mod_Generation_Fruits_Step: fitness ~ Generation + Step + Fruit_s + (1 | Generation)
## mod_Double_gene: fitness ~ Generation * Fruit_s + (1 | Generation)
## mod_Doubleinteraction: fitness ~ Generation + Step + Fruit_s + Fruit_s:Generation + 
## mod_Doubleinteraction:     (1 | Generation)
## mod_Double: fitness ~ Step * Fruit_s + (1 | Generation)
## mod_Double_Generation: fitness ~ Generation + Step * Fruit_s + (1 | Generation)
## mod_Tripleinteraction: fitness ~ Generation * Step * Fruit_s + (1 | Generation)
## mod_Lines: fitness ~ Line + (1 | Generation)
## mod_DoublePopulationGeneration: fitness ~ Generation * Line + (1 | Generation)
##                                npar    AIC    BIC  logLik deviance   Chisq Df
## mod_nul                           3 336.65 347.18 -165.33   330.65           
## mod_Generation                    4 323.01 337.05 -157.50   315.01 15.6458  1
## mod_Fruits                        5 321.74 339.29 -155.87   311.74  3.2688  1
## mod_Step                          5 305.39 322.93 -147.69   295.39 16.3510  0
## mod_Generation_Step               6 305.95 327.01 -146.98   293.95  1.4362  1
## mod_Generation_Fruits             6 308.84 329.90 -148.42   296.84  0.0000  0
## mod_Fruits_Step                   7 291.48 316.05 -138.74   277.48 19.3583  1
## mod_Generation_Fruits_Step        8 291.95 320.02 -137.97   275.95  1.5325  1
## mod_Double_gene                   8 306.40 334.48 -145.20   290.40  0.0000  0
## mod_Doubleinteraction            10 289.56 324.65 -134.78   269.56 20.8451  2
## mod_Double                       11 291.19 329.79 -134.59   269.19  0.3679  1
## mod_Double_Generation            12 291.62 333.73 -133.81   267.62  1.5683  1
## mod_Tripleinteraction            20 290.03 360.22 -125.02   250.03 17.5887  8
## mod_Lines                        31 304.85 413.64 -121.42   242.85  7.1828 11
## mod_DoublePopulationGeneration   60 322.24 532.80 -101.12   202.24 40.6128 29
##                                Pr(>Chisq)    
## mod_nul                                      
## mod_Generation                  7.638e-05 ***
## mod_Fruits                        0.07061 .  
## mod_Step                        < 2.2e-16 ***
## mod_Generation_Step               0.23076    
## mod_Generation_Fruits             1.00000    
## mod_Fruits_Step                 1.083e-05 ***
## mod_Generation_Fruits_Step        0.21574    
## mod_Double_gene                   1.00000    
## mod_Doubleinteraction           2.975e-05 ***
## mod_Double                        0.54414    
## mod_Double_Generation             0.21045    
## mod_Tripleinteraction             0.02453 *  
## mod_Lines                         0.78410    
## mod_DoublePopulationGeneration    0.07441 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod_Fruits_Step)
## Linear mixed model fit by REML ['lmerMod']
## Formula: fitness ~ Step + Fruit_s + (1 | Generation)
##    Data: data_sum
## Weights: N
## 
## REML criterion at convergence: 294.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8621 -0.4273  0.0855  0.5085  2.5738 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Generation (Intercept) 0.03646  0.1909  
##  Residual               3.14218  1.7726  
## Number of obs: 247, groups:  Generation, 24
## 
## Fixed effects:
##                     Estimate Std. Error t value
## (Intercept)         -0.33639    0.12545  -2.681
## Steppool             0.55900    0.16167   3.458
## Stepsecond_postpool  1.00604    0.13160   7.645
## Fruit_sCranberry    -0.15931    0.05744  -2.774
## Fruit_sStrawberry   -0.26942    0.06258  -4.305
## 
## Correlation of Fixed Effects:
##             (Intr) Steppl Stpsc_ Frt_sC
## Steppool    -0.689                     
## Stpscnd_pst -0.864  0.670              
## Frt_sCrnbrr -0.278 -0.043 -0.005       
## Frt_sStrwbr -0.271 -0.036  0.013  0.591
#Posthoc
mod_Best <- lme4::lmer(fitness ~ Step -1 + Fruit_s + (1|Generation),
                     weights = N, data = data_sum)
emmeans::emmeans(mod_Best, list(pairwise ~ Step:Fruit_s), adjust = "tukey") #
## $`emmeans of Step, Fruit_s`
##  Step            Fruit_s     emmean     SE    df lower.CL upper.CL
##  first_prepool   Cherry     -0.3364 0.1285   300   -0.694   0.0216
##  pool            Cherry      0.2226 0.1181 83203   -0.104   0.5493
##  second_postpool Cherry      0.6697 0.0678  2115    0.482   0.8574
##  first_prepool   Cranberry  -0.4957 0.1256   287   -0.846  -0.1456
##  pool            Cranberry   0.0633 0.1115 92051   -0.245   0.3718
##  second_postpool Cranberry   0.5103 0.0618  1382    0.339   0.6814
##  first_prepool   Strawberry -0.6058 0.1271   294   -0.960  -0.2516
##  pool            Strawberry -0.0468 0.1136 91454   -0.361   0.2672
##  second_postpool Strawberry  0.4002 0.0669  1974    0.215   0.5854
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 9 estimates 
## 
## $`pairwise differences of Step, Fruit_s`
##  contrast                                               estimate     SE     df
##  first_prepool Cherry - pool Cherry                       -0.559 0.1640    797
##  first_prepool Cherry - second_postpool Cherry            -1.006 0.1347    303
##  first_prepool Cherry - first_prepool Cranberry            0.159 0.0575 114310
##  first_prepool Cherry - pool Cranberry                    -0.400 0.1715    920
##  first_prepool Cherry - second_postpool Cranberry         -0.847 0.1463    407
##  first_prepool Cherry - first_prepool Strawberry           0.269 0.0626 112175
##  first_prepool Cherry - pool Strawberry                   -0.290 0.1735    987
##  first_prepool Cherry - second_postpool Strawberry        -0.737 0.1493    455
##  pool Cherry - second_postpool Cherry                     -0.447 0.1227  12246
##  pool Cherry - first_prepool Cranberry                     0.718 0.1761   1109
##  pool Cherry - pool Cranberry                              0.159 0.0575 114310
##  pool Cherry - second_postpool Cranberry                  -0.288 0.1382  17293
##  pool Cherry - first_prepool Strawberry                    0.828 0.1776   1124
##  pool Cherry - pool Strawberry                             0.269 0.0626 112175
##  pool Cherry - second_postpool Strawberry                 -0.178 0.1411  19292
##  second_postpool Cherry - first_prepool Cranberry          1.165 0.1467    454
##  second_postpool Cherry - pool Cranberry                   0.606 0.1327  17931
##  second_postpool Cherry - second_postpool Cranberry        0.159 0.0575 114310
##  second_postpool Cherry - first_prepool Strawberry         1.275 0.1478    457
##  second_postpool Cherry - pool Strawberry                  0.716 0.1343  18072
##  second_postpool Cherry - second_postpool Strawberry       0.269 0.0626 112175
##  first_prepool Cranberry - pool Cranberry                 -0.559 0.1640    797
##  first_prepool Cranberry - second_postpool Cranberry      -1.006 0.1347    303
##  first_prepool Cranberry - first_prepool Strawberry        0.110 0.0545 141545
##  first_prepool Cranberry - pool Strawberry                -0.449 0.1731   1016
##  first_prepool Cranberry - second_postpool Strawberry     -0.896 0.1463    439
##  pool Cranberry - second_postpool Cranberry               -0.447 0.1227  12246
##  pool Cranberry - first_prepool Strawberry                 0.669 0.1726    962
##  pool Cranberry - pool Strawberry                          0.110 0.0545 141545
##  pool Cranberry - second_postpool Strawberry              -0.337 0.1350  18114
##  second_postpool Cranberry - first_prepool Strawberry      1.116 0.1444    394
##  second_postpool Cranberry - pool Strawberry               0.557 0.1335  16219
##  second_postpool Cranberry - second_postpool Strawberry    0.110 0.0545 141545
##  first_prepool Strawberry - pool Strawberry               -0.559 0.1640    797
##  first_prepool Strawberry - second_postpool Strawberry    -1.006 0.1347    303
##  pool Strawberry - second_postpool Strawberry             -0.447 0.1227  12246
##  t.ratio p.value
##  -3.408  0.0197 
##  -7.467  <.0001 
##   2.773  0.1235 
##  -2.330  0.3250 
##  -5.789  <.0001 
##   4.303  0.0006 
##  -1.669  0.7656 
##  -4.933  <.0001 
##  -3.644  0.0083 
##   4.080  0.0016 
##   2.773  0.1235 
##  -2.083  0.4853 
##   4.664  0.0001 
##   4.303  0.0006 
##  -1.259  0.9430 
##   7.946  <.0001 
##   4.568  0.0002 
##   2.773  0.1235 
##   8.629  <.0001 
##   5.335  <.0001 
##   4.303  0.0006 
##  -3.408  0.0197 
##  -7.467  <.0001 
##   2.020  0.5293 
##  -2.593  0.1905 
##  -6.124  <.0001 
##  -3.644  0.0083 
##   3.876  0.0036 
##   2.020  0.5293 
##  -2.496  0.2340 
##   7.731  <.0001 
##   4.174  0.0010 
##   2.020  0.5293 
##  -3.408  0.0197 
##  -7.467  <.0001 
##  -3.644  0.0083 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 9 estimates
emmeans::emmeans(mod_Best, list(pairwise ~ Fruit_s), adjust = "tukey") #
## $`emmeans of Fruit_s`
##  Fruit_s     emmean     SE   df lower.CL upper.CL
##  Cherry      0.1853 0.0708 1709   0.0161   0.3545
##  Cranberry   0.0260 0.0632 1231  -0.1252   0.1772
##  Strawberry -0.0841 0.0671 1475  -0.2444   0.0762
## 
## Results are averaged over the levels of: Step 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Fruit_s`
##  contrast               estimate     SE     df t.ratio p.value
##  Cherry - Cranberry        0.159 0.0575 114310 2.773   0.0154 
##  Cherry - Strawberry       0.269 0.0626 112175 4.303   <.0001 
##  Cranberry - Strawberry    0.110 0.0545 141545 2.020   0.1074 
## 
## Results are averaged over the levels of: Step 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates
emmeans::emmeans(mod_Best, list(pairwise ~ Step), adjust = "tukey") #
## $`emmeans of Step`
##  Step             emmean     SE    df lower.CL upper.CL
##  first_prepool   -0.4793 0.1226   252   -0.774   -0.185
##  pool             0.0797 0.1094 81509   -0.181    0.341
##  second_postpool  0.5267 0.0562   970    0.392    0.661
## 
## Results are averaged over the levels of: Fruit_s 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Step`
##  contrast                        estimate    SE    df t.ratio p.value
##  first_prepool - pool              -0.559 0.164   797 -3.408  0.0020 
##  first_prepool - second_postpool   -1.006 0.135   303 -7.467  <.0001 
##  pool - second_postpool            -0.447 0.123 12246 -3.644  0.0008 
## 
## Results are averaged over the levels of: Fruit_s 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates
summary(mod_Best)
## Linear mixed model fit by REML ['lmerMod']
## Formula: fitness ~ Step - 1 + Fruit_s + (1 | Generation)
##    Data: data_sum
## Weights: N
## 
## REML criterion at convergence: 294.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8621 -0.4273  0.0855  0.5085  2.5738 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Generation (Intercept) 0.03646  0.1909  
##  Residual               3.14218  1.7726  
## Number of obs: 247, groups:  Generation, 24
## 
## Fixed effects:
##                     Estimate Std. Error t value
## Stepfirst_prepool   -0.33639    0.12545  -2.681
## Steppool             0.22261    0.11808   1.885
## Stepsecond_postpool  0.66965    0.06728   9.953
## Fruit_sCranberry    -0.15931    0.05744  -2.774
## Fruit_sStrawberry   -0.26942    0.06258  -4.305
## 
## Correlation of Fixed Effects:
##             Stpfr_ Steppl Stpsc_ Frt_sC
## Steppool     0.120                     
## Stpscnd_pst  0.175  0.220              
## Frt_sCrnbrr -0.278 -0.355 -0.528       
## Frt_sStrwbr -0.271 -0.337 -0.480  0.591
################### BELOW TEMP: A SUPR
############### Percent of Increase in fitness
Calcul_increase_fitness <- expand.grid(Step = unique(data_sum$Step),
                                     Fruit_s = unique(data_sum$Fruit_s))

Calcul_increase_fitness$fitness_predicted <- predict(mod_Best, newdata = Calcul_increase_fitness, 
                          re.form= NA, type = "response")

((Calcul_increase_fitness$fitness_predicted[Calcul_increase_fitness$Step=="second_postpool"]-
    Calcul_increase_fitness$fitness_predicted[Calcul_increase_fitness$Step=="first_prepool"])/
  abs(Calcul_increase_fitness$fitness_predicted[Calcul_increase_fitness$Step=="first_prepool"]))*100
## [1] 299.0683 202.9538 166.0655
Calcul_increase_fitness
##              Step    Fruit_s fitness_predicted
## 1   first_prepool     Cherry       -0.33639290
## 2 second_postpool     Cherry        0.66965164
## 3            pool     Cherry        0.22261164
## 4   first_prepool  Cranberry       -0.49570136
## 5 second_postpool  Cranberry        0.51034318
## 6            pool  Cranberry        0.06330318
## 7   first_prepool Strawberry       -0.60581186
## 8 second_postpool Strawberry        0.40023267
## 9            pool Strawberry       -0.04680732

2.2 Plot

#Position dodge
pd <- ggplot2::position_dodge(0.3) # move them .05 to the left and right

#Extract slope and intercept
dat_predict_allfruits <- expand.grid(Generation =as.numeric(levels(as.factor(data_sum$Generation))),
                                     Fruit_s=unique(data_sum$Fruit_s))

dat_predict_allfruits$Step <- ifelse(dat_predict_allfruits$Generation<= 7, "first_prepool",ifelse(dat_predict_allfruits$Generation>=12, "second_postpool", "pool"))

dat_predict_allfruits$fitness_predicted <- predict(mod_Best, newdata = dat_predict_allfruits, 
                          re.form= NA, type = "response")

#REAL DATA
#Add G1 G6 and G7
TEMP_lineG1 <- c(rep(NA, 2), 1,rep(NA, 5))
TEMP_lineG6 <- c(rep(NA, 2), 6,rep(NA, 5))
TEMP_lineG7 <- c(rep(NA, 2), 7,rep(NA, 5))

TEMP_total <- rbind(data_sum,TEMP_lineG1,TEMP_lineG1,TEMP_lineG1,TEMP_lineG6,TEMP_lineG6,TEMP_lineG6,TEMP_lineG7,TEMP_lineG7,TEMP_lineG7)
TEMP_total$Fruit_s[TEMP_total$Generation == "6"] <- c("Strawberry", "Cranberry", "Cherry")
TEMP_total$Fruit_s[TEMP_total$Generation == "1"] <- c("Strawberry", "Cranberry", "Cherry")
TEMP_total$Fruit_s[TEMP_total$Generation == "7"] <- c("Strawberry", "Cranberry", "Cherry")
tail(TEMP_total)
##     Line    Fruit_s Generation Step  N Nb_adults sd fitness se_fitness
## 251 <NA> Strawberry          6 <NA> NA        NA NA      NA         NA
## 252 <NA>  Cranberry          6 <NA> NA        NA NA      NA         NA
## 253 <NA>     Cherry          6 <NA> NA        NA NA      NA         NA
## 254 <NA> Strawberry          7 <NA> NA        NA NA      NA         NA
## 255 <NA>  Cranberry          7 <NA> NA        NA NA      NA         NA
## 256 <NA>     Cherry          7 <NA> NA        NA NA      NA         NA
## Add label 
TEMP_anno <- data.frame(x1 = c(3.5, 3.5, 10, 3.5, 3.5, 10, 3.5, 3.5, 10), 
                        x2 = c(9, 17.5, 17.5, 9, 17.5, 17.5, 9, 17.5, 17.5),
                        y1 = c(1, 2, 1.2, 1, 2, 1.2, 1, 2, 1.2), 
                        y2 = c(1.25, 2.25, 1.45, 1.25, 2.25, 1.45, 1.25, 2.25, 1.45),
                        xstar = c(6.5, 10, 14, 6.5, 10, 14, 6.5, 10, 14), 
                        ystar = c(1.5, 2.5, 1.7, 1.5, 2.5, 1.7, 1.5, 2.5, 1.7),
                        lab = c("**", "***", "***",  "**", "***", "***",  "**", "***", "***"),
                        Fruit_s = c("Cherry", "Cherry", "Cherry",
                               "Cranberry", "Cranberry", "Cranberry",
                               "Strawberry", "Strawberry", "Strawberry"), 
                        Line = NA)
TEMP_anno
##     x1   x2  y1   y2 xstar ystar lab    Fruit_s Line
## 1  3.5  9.0 1.0 1.25   6.5   1.5  **     Cherry   NA
## 2  3.5 17.5 2.0 2.25  10.0   2.5 ***     Cherry   NA
## 3 10.0 17.5 1.2 1.45  14.0   1.7 ***     Cherry   NA
## 4  3.5  9.0 1.0 1.25   6.5   1.5  **  Cranberry   NA
## 5  3.5 17.5 2.0 2.25  10.0   2.5 ***  Cranberry   NA
## 6 10.0 17.5 1.2 1.45  14.0   1.7 ***  Cranberry   NA
## 7  3.5  9.0 1.0 1.25   6.5   1.5  ** Strawberry   NA
## 8  3.5 17.5 2.0 2.25  10.0   2.5 *** Strawberry   NA
## 9 10.0 17.5 1.2 1.45  14.0   1.7 *** Strawberry   NA
TEMP_title <- data.frame(xtitle = c(14, 14, 14),
                       ytitle = c(3, 3, 3),
                       title = c("Cherry", "Cranberry", "Strawberry"),
                       Fruit_s = c("Cherry", "Cranberry", "Strawberry"), 
                       Line = NA)
TEMP_title
##   xtitle ytitle      title    Fruit_s Line
## 1     14      3     Cherry     Cherry   NA
## 2     14      3  Cranberry  Cranberry   NA
## 3     14      3 Strawberry Strawberry   NA
PLOT_FITNESS_CHERRY <- ggplot(data = TEMP_total[TEMP_total$Fruit_s == "Cherry",], 
                            aes(x = factor(Generation),group = Line, y = fitness, colour =Fruit_s)) + 
  geom_errorbar(aes(ymin =fitness-1.96*se_fitness, ymax = fitness + 1.96*se_fitness),
                width=.1,position = pd, size = 0.2,color = "black") + 
  geom_line(size = 0.3,position = pd) + 
  geom_line(data = dat_predict_allfruits[dat_predict_allfruits$Fruit_s == "Cherry",], 
                     aes(x = factor(Generation), y = fitness_predicted,
                         colour = "black", group = Step), size = 0.5) + 
  geom_point(size =1, position = pd, shape =21, fill = "white") + 
  ylim(-3, 3.05) + 
  ylab("Fitness") + 
  xlab("Generation") + 
  geom_text(data = TEMP_anno[TEMP_anno$Fruit_s == "Cherry",], aes(x = xstar,  y = ystar, label = lab), size =3.3) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cherry",], aes(x = x1, xend = x1, 
           y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cherry",], aes(x = x2, xend = x2, 
           y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cherry",], aes(x = x1, xend = x2, 
           y = y2, yend = y2), size = 0.4) + 
  scale_color_manual(values = c("black", "#BC3C6D", "#FDB424", "#3FAA96")) + 
  ggtitle("Cherry") + 
  theme_LO_adaptation + theme(plot.title = element_text(color = "#BC3C6D"))
PLOT_FITNESS_CHERRY

PLOT_FITNESS_CRANB <- ggplot2::ggplot(data = TEMP_total[TEMP_total$Fruit_s == "Cranberry",], 
                            aes(x = factor(Generation),group = Line, y = fitness, colour =Fruit_s)) + 
  geom_errorbar(aes(ymin =fitness-1.96*se_fitness, ymax = fitness + 1.96*se_fitness),
                width=.1,position = pd, size = 0.2,color = "black") + 
  geom_line(size = 0.3,position = pd) + 
  geom_line(data = dat_predict_allfruits[dat_predict_allfruits$Fruit_s == "Cranberry",], aes(x = factor(Generation), y = fitness_predicted, colour = "black", group = Step),
                size = 0.5) + 
  geom_point(size =1, position = pd, shape =21, fill = "white") + 
  ylim(-3, 3.05) + 
  ylab("Fitness") + 
  xlab("Generation") + 
  geom_text(data = TEMP_anno[TEMP_anno$Fruit_s == "Cranberry",], aes(x = xstar,  y = ystar, label = lab), size =3.3) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cranberry",], aes(x = x1, xend = x1, 
           y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cranberry",], aes(x = x2, xend = x2, 
           y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cranberry",], aes(x = x1, xend = x2, 
           y = y2, yend = y2), size = 0.4) + 
  scale_color_manual(values = c("black", "#FDB424")) + 
  ggtitle("Cranberry") + 
  theme_LO_adaptation + theme(plot.title = element_text(color = "#FDB424"))
PLOT_FITNESS_CRANB

PLOT_FITNESS_STRAW <- ggplot2::ggplot(data = TEMP_total[TEMP_total$Fruit_s == "Strawberry",], 
                            aes(x = factor(Generation),group = Line, y = fitness, colour =Fruit_s)) + 
  geom_errorbar(aes(ymin =fitness-1.96*se_fitness, ymax = fitness + 1.96*se_fitness),
                width=.1,position = pd, size = 0.2,color = "black") + 
  geom_line(size = 0.3,position = pd) + 
  geom_line(data = dat_predict_allfruits[dat_predict_allfruits$Fruit_s == "Strawberry",], 
            aes(x = factor(Generation), y = fitness_predicted, colour = "black", group = Step),
            size = 0.5) + 
  geom_point(size =1, position = pd, shape =21, fill = "white") + 
  ylim(-3, 3.05) + 
  ylab("Fitness") + 
  xlab("Generation") + 
  geom_text(data = TEMP_anno[TEMP_anno$Fruit_s == "Strawberry",], 
            aes(x = xstar,  y = ystar, label = lab), size =3.3) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Strawberry",], 
               aes(x = x1, xend = x1, y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Strawberry",], 
               aes(x = x2, xend = x2, y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Strawberry",], 
               aes(x = x1, xend = x2, y = y2, yend = y2), size = 0.4) + 
  scale_color_manual(values = c("black", "#3FAA96")) + 
  ggtitle("Strawberry") + 
  theme_LO_adaptation + theme(plot.title = element_text(color = "#3FAA96"))
PLOT_FITNESS_STRAW

DYNAMIQUE_THREE_JOIN <- cowplot::ggdraw() + 
  cowplot::draw_plot(PLOT_FITNESS_CHERRY + theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank()), 
            x = 0, y = 0.66, width = 1, height = 0.33) + 
  cowplot::draw_plot(PLOT_FITNESS_CRANB + theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank()), 
            x = 0, y = 0.33, width = 1, height = 0.33) + 
  cowplot::draw_plot(PLOT_FITNESS_STRAW, 
            x = 0, y = 00, width = 1, height = 0.33) + 
  cowplot::draw_plot_label(c("A", "B", "C"),  
                  x = c(0, 0, 0), 
                  y = c(1, 0.66, 0.33), 
                  hjust = c(-0.5, -0.5, -0.5), 
                  vjust = c(1.5, 1.5, 1.5),
                  size = 12) 
 DYNAMIQUE_THREE_JOIN

cowplot::save_plot(file =here::here("figures", "FIG_Adaptation.pdf"), DYNAMIQUE_THREE_JOIN, base_height = 17/cm(1), base_width = 11/cm(1), dpi = 1200)

2.3 Analyses with indic

#Posthoc
mod_Best <- lme4::lmer(fitness ~ Step -1 + Fruit_s + (1|Generation),
                     weights = N, data = data_sum)
emmeans::emmeans(mod_Best, list(pairwise ~ Step:Fruit_s), adjust = "tukey") #
## $`emmeans of Step, Fruit_s`
##  Step            Fruit_s     emmean     SE    df lower.CL upper.CL
##  first_prepool   Cherry     -0.3364 0.1285   300   -0.694   0.0216
##  pool            Cherry      0.2226 0.1181 83203   -0.104   0.5493
##  second_postpool Cherry      0.6697 0.0678  2115    0.482   0.8574
##  first_prepool   Cranberry  -0.4957 0.1256   287   -0.846  -0.1456
##  pool            Cranberry   0.0633 0.1115 92051   -0.245   0.3718
##  second_postpool Cranberry   0.5103 0.0618  1382    0.339   0.6814
##  first_prepool   Strawberry -0.6058 0.1271   294   -0.960  -0.2516
##  pool            Strawberry -0.0468 0.1136 91454   -0.361   0.2672
##  second_postpool Strawberry  0.4002 0.0669  1974    0.215   0.5854
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 9 estimates 
## 
## $`pairwise differences of Step, Fruit_s`
##  contrast                                               estimate     SE     df
##  first_prepool Cherry - pool Cherry                       -0.559 0.1640    797
##  first_prepool Cherry - second_postpool Cherry            -1.006 0.1347    303
##  first_prepool Cherry - first_prepool Cranberry            0.159 0.0575 114310
##  first_prepool Cherry - pool Cranberry                    -0.400 0.1715    920
##  first_prepool Cherry - second_postpool Cranberry         -0.847 0.1463    407
##  first_prepool Cherry - first_prepool Strawberry           0.269 0.0626 112175
##  first_prepool Cherry - pool Strawberry                   -0.290 0.1735    987
##  first_prepool Cherry - second_postpool Strawberry        -0.737 0.1493    455
##  pool Cherry - second_postpool Cherry                     -0.447 0.1227  12246
##  pool Cherry - first_prepool Cranberry                     0.718 0.1761   1109
##  pool Cherry - pool Cranberry                              0.159 0.0575 114310
##  pool Cherry - second_postpool Cranberry                  -0.288 0.1382  17293
##  pool Cherry - first_prepool Strawberry                    0.828 0.1776   1124
##  pool Cherry - pool Strawberry                             0.269 0.0626 112175
##  pool Cherry - second_postpool Strawberry                 -0.178 0.1411  19292
##  second_postpool Cherry - first_prepool Cranberry          1.165 0.1467    454
##  second_postpool Cherry - pool Cranberry                   0.606 0.1327  17931
##  second_postpool Cherry - second_postpool Cranberry        0.159 0.0575 114310
##  second_postpool Cherry - first_prepool Strawberry         1.275 0.1478    457
##  second_postpool Cherry - pool Strawberry                  0.716 0.1343  18072
##  second_postpool Cherry - second_postpool Strawberry       0.269 0.0626 112175
##  first_prepool Cranberry - pool Cranberry                 -0.559 0.1640    797
##  first_prepool Cranberry - second_postpool Cranberry      -1.006 0.1347    303
##  first_prepool Cranberry - first_prepool Strawberry        0.110 0.0545 141545
##  first_prepool Cranberry - pool Strawberry                -0.449 0.1731   1016
##  first_prepool Cranberry - second_postpool Strawberry     -0.896 0.1463    439
##  pool Cranberry - second_postpool Cranberry               -0.447 0.1227  12246
##  pool Cranberry - first_prepool Strawberry                 0.669 0.1726    962
##  pool Cranberry - pool Strawberry                          0.110 0.0545 141545
##  pool Cranberry - second_postpool Strawberry              -0.337 0.1350  18114
##  second_postpool Cranberry - first_prepool Strawberry      1.116 0.1444    394
##  second_postpool Cranberry - pool Strawberry               0.557 0.1335  16219
##  second_postpool Cranberry - second_postpool Strawberry    0.110 0.0545 141545
##  first_prepool Strawberry - pool Strawberry               -0.559 0.1640    797
##  first_prepool Strawberry - second_postpool Strawberry    -1.006 0.1347    303
##  pool Strawberry - second_postpool Strawberry             -0.447 0.1227  12246
##  t.ratio p.value
##  -3.408  0.0197 
##  -7.467  <.0001 
##   2.773  0.1235 
##  -2.330  0.3250 
##  -5.789  <.0001 
##   4.303  0.0006 
##  -1.669  0.7656 
##  -4.933  <.0001 
##  -3.644  0.0083 
##   4.080  0.0016 
##   2.773  0.1235 
##  -2.083  0.4853 
##   4.664  0.0001 
##   4.303  0.0006 
##  -1.259  0.9430 
##   7.946  <.0001 
##   4.568  0.0002 
##   2.773  0.1235 
##   8.629  <.0001 
##   5.335  <.0001 
##   4.303  0.0006 
##  -3.408  0.0197 
##  -7.467  <.0001 
##   2.020  0.5293 
##  -2.593  0.1905 
##  -6.124  <.0001 
##  -3.644  0.0083 
##   3.876  0.0036 
##   2.020  0.5293 
##  -2.496  0.2340 
##   7.731  <.0001 
##   4.174  0.0010 
##   2.020  0.5293 
##  -3.408  0.0197 
##  -7.467  <.0001 
##  -3.644  0.0083 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 9 estimates
emmeans::emmeans(mod_Best, list(pairwise ~ Fruit_s), adjust = "tukey") #
## $`emmeans of Fruit_s`
##  Fruit_s     emmean     SE   df lower.CL upper.CL
##  Cherry      0.1853 0.0708 1709   0.0161   0.3545
##  Cranberry   0.0260 0.0632 1231  -0.1252   0.1772
##  Strawberry -0.0841 0.0671 1475  -0.2444   0.0762
## 
## Results are averaged over the levels of: Step 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Fruit_s`
##  contrast               estimate     SE     df t.ratio p.value
##  Cherry - Cranberry        0.159 0.0575 114310 2.773   0.0154 
##  Cherry - Strawberry       0.269 0.0626 112175 4.303   <.0001 
##  Cranberry - Strawberry    0.110 0.0545 141545 2.020   0.1074 
## 
## Results are averaged over the levels of: Step 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates
emmeans::emmeans(mod_Best, list(pairwise ~ Step), adjust = "tukey") #
## $`emmeans of Step`
##  Step             emmean     SE    df lower.CL upper.CL
##  first_prepool   -0.4793 0.1226   252   -0.774   -0.185
##  pool             0.0797 0.1094 81509   -0.181    0.341
##  second_postpool  0.5267 0.0562   970    0.392    0.661
## 
## Results are averaged over the levels of: Fruit_s 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Step`
##  contrast                        estimate    SE    df t.ratio p.value
##  first_prepool - pool              -0.559 0.164   797 -3.408  0.0020 
##  first_prepool - second_postpool   -1.006 0.135   303 -7.467  <.0001 
##  pool - second_postpool            -0.447 0.123 12246 -3.644  0.0008 
## 
## Results are averaged over the levels of: Fruit_s 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates
summary(mod_Best)
## Linear mixed model fit by REML ['lmerMod']
## Formula: fitness ~ Step - 1 + Fruit_s + (1 | Generation)
##    Data: data_sum
## Weights: N
## 
## REML criterion at convergence: 294.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8621 -0.4273  0.0855  0.5085  2.5738 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Generation (Intercept) 0.03646  0.1909  
##  Residual               3.14218  1.7726  
## Number of obs: 247, groups:  Generation, 24
## 
## Fixed effects:
##                     Estimate Std. Error t value
## Stepfirst_prepool   -0.33639    0.12545  -2.681
## Steppool             0.22261    0.11808   1.885
## Stepsecond_postpool  0.66965    0.06728   9.953
## Fruit_sCranberry    -0.15931    0.05744  -2.774
## Fruit_sStrawberry   -0.26942    0.06258  -4.305
## 
## Correlation of Fixed Effects:
##             Stpfr_ Steppl Stpsc_ Frt_sC
## Steppool     0.120                     
## Stpscnd_pst  0.175  0.220              
## Frt_sCrnbrr -0.278 -0.355 -0.528       
## Frt_sStrwbr -0.271 -0.337 -0.480  0.591
############## 0 Add indic 
data_sum$indic_phase1 <- ifelse(data_sum$Step=="first_prepool", 1, 0)
data_sum$indic_phase2 <- ifelse(data_sum$Step=="pool", 1, 0)
data_sum$indic_phase3 <- ifelse(data_sum$Step=="second_postpool", 1, 0)
data_sum$indic_cherry<- ifelse(data_sum$Fruit_s=="Cherry", "Cherry", "Others")
data_sum$indic_cranberry <- ifelse(data_sum$Fruit_s=="Cranberry", "Cranberry", "Others")
data_sum$indic_strawberry <- ifelse(data_sum$Fruit_s=="Strawberry",  "Strawberry", "Others")


############## 1 Test indic sur step * Fruit 
mod_Fruit_Indic1 <- lme4::lmer(fitness ~ Fruit_s*indic_phase1 + (1|Generation),
            weights = N, data = data_sum)

mod_Fruit_Indic2 <- lme4::lmer(fitness ~ Fruit_s*indic_phase2 + (1|Generation),
            weights = N, data = data_sum)

mod_Fruit_Indic3 <- lme4::lmer(fitness ~ Fruit_s*indic_phase3 + (1|Generation),
            weights = N, data = data_sum)


############## 2 Test indic sur fRuit * Fruit 
mod_Cherry_Step <- lme4::lmer(fitness ~ Step*indic_cherry + (1|Generation),
            weights = N, data = data_sum)

mod_Cranberry_Step <- lme4::lmer(fitness ~ Step*indic_cranberry + (1|Generation),
            weights = N, data = data_sum)

mod_Strawberry_Step <- lme4::lmer(fitness ~ Step*indic_strawberry + (1|Generation),
            weights = N, data = data_sum)


############## 3 Test indic sur fRuit * Test indic phase 3
mod_Cherry_indic1 <- lme4::lmer(fitness ~ indic_phase1*indic_cherry + (1|Generation),
            weights = N, data = data_sum)
mod_Cherry_indic2 <- lme4::lmer(fitness ~ indic_phase2*indic_cherry + (1|Generation),
            weights = N, data = data_sum)
mod_Cherry_indic3 <- lme4::lmer(fitness ~ indic_phase3*indic_cherry + (1|Generation),
            weights = N, data = data_sum)


MuMIn::AICc(mod_Fruit_Indic1,mod_Fruit_Indic2,mod_Fruit_Indic3,
            mod_Cherry_Step,mod_Cranberry_Step,mod_Strawberry_Step,
            mod_Cherry_indic1,mod_Cherry_indic2,mod_Cherry_indic3)
##                     df     AICc
## mod_Fruit_Indic1     8 314.8545
## mod_Fruit_Indic2     8 342.1753
## mod_Fruit_Indic3     8 321.9248
## mod_Cherry_Step      8 307.8003
## mod_Cranberry_Step   8 328.4648
## mod_Strawberry_Step  8 316.4444
## mod_Cherry_indic1    6 309.9783
## mod_Cherry_indic2    6 336.5756
## mod_Cherry_indic3    6 315.4998
anova(mod_Fruit_Indic1,mod_Fruit_Indic2,mod_Fruit_Indic3,
            mod_Cherry_Step,mod_Cranberry_Step,mod_Strawberry_Step,
            mod_Cherry_indic1,mod_Cherry_indic2,mod_Cherry_indic3)
## Data: data_sum
## Models:
## mod_Cherry_indic1: fitness ~ indic_phase1 * indic_cherry + (1 | Generation)
## mod_Cherry_indic2: fitness ~ indic_phase2 * indic_cherry + (1 | Generation)
## mod_Cherry_indic3: fitness ~ indic_phase3 * indic_cherry + (1 | Generation)
## mod_Fruit_Indic1: fitness ~ Fruit_s * indic_phase1 + (1 | Generation)
## mod_Fruit_Indic2: fitness ~ Fruit_s * indic_phase2 + (1 | Generation)
## mod_Fruit_Indic3: fitness ~ Fruit_s * indic_phase3 + (1 | Generation)
## mod_Cherry_Step: fitness ~ Step * indic_cherry + (1 | Generation)
## mod_Cranberry_Step: fitness ~ Step * indic_cranberry + (1 | Generation)
## mod_Strawberry_Step: fitness ~ Step * indic_strawberry + (1 | Generation)
##                     npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)
## mod_Cherry_indic1      6 297.91 318.96 -142.95   285.91                      
## mod_Cherry_indic2      6 326.22 347.28 -157.11   314.22  0.0000  0   1.000000
## mod_Cherry_indic3      6 302.24 323.29 -145.12   290.24 23.9844  0  < 2.2e-16
## mod_Fruit_Indic1       8 296.71 324.79 -140.36   280.71  9.5266  2   0.008537
## mod_Fruit_Indic2       8 325.31 353.38 -154.65   309.31  0.0000  0   1.000000
## mod_Fruit_Indic3       8 301.90 329.98 -142.95   285.90 23.4071  0  < 2.2e-16
## mod_Cherry_Step        8 290.29 318.37 -137.15   274.29 11.6067  0  < 2.2e-16
## mod_Cranberry_Step     8 310.11 338.19 -147.06   294.11  0.0000  0   1.000000
## mod_Strawberry_Step    8 298.37 326.45 -141.19   282.37 11.7394  0  < 2.2e-16
##                        
## mod_Cherry_indic1      
## mod_Cherry_indic2      
## mod_Cherry_indic3   ***
## mod_Fruit_Indic1    ** 
## mod_Fruit_Indic2       
## mod_Fruit_Indic3    ***
## mod_Cherry_Step     ***
## mod_Cranberry_Step     
## mod_Strawberry_Step ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Posthoc
mod_Best <- lme4::lmer(fitness ~ Step*indic_cherry + (1|Generation),
            weights = N, data = data_sum)

emmeans::emmeans(mod_Best, list(pairwise ~ Step:indic_cherry), adjust = "tukey") #
## $`emmeans of Step, indic_cherry`
##  Step            indic_cherry    emmean     SE     df lower.CL upper.CL
##  first_prepool   Cherry       -0.607592 0.1631    674   -1.038   -0.177
##  pool            Cherry        0.321562 0.1674 111144   -0.119    0.762
##  second_postpool Cherry        0.701026 0.0712   2598    0.514    0.889
##  first_prepool   Others       -0.442650 0.1288    340   -0.783   -0.102
##  pool            Others        0.000256 0.1106 107399   -0.291    0.291
##  second_postpool Others        0.455445 0.0585   1121    0.301    0.610
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 6 estimates 
## 
## $`pairwise differences of Step, indic_cherry`
##  contrast                                        estimate     SE     df t.ratio
##  first_prepool Cherry - pool Cherry                -0.929 0.2337   2618 -3.976 
##  first_prepool Cherry - second_postpool Cherry     -1.309 0.1779    796 -7.354 
##  first_prepool Cherry - first_prepool Others       -0.165 0.1484  27353 -1.111 
##  first_prepool Cherry - pool Others                -0.608 0.1970   1380 -3.085 
##  first_prepool Cherry - second_postpool Others     -1.063 0.1733    710 -6.135 
##  pool Cherry - second_postpool Cherry              -0.379 0.1819  38241 -2.086 
##  pool Cherry - first_prepool Others                 0.764 0.2112   2298  3.619 
##  pool Cherry - pool Others                          0.321 0.1479 504858  2.173 
##  pool Cherry - second_postpool Others              -0.134 0.1773  34407 -0.755 
##  second_postpool Cherry - first_prepool Others      1.144 0.1471    475  7.773 
##  second_postpool Cherry - pool Others               0.701 0.1315  17346  5.328 
##  second_postpool Cherry - second_postpool Others    0.246 0.0614 125672  3.997 
##  first_prepool Others - pool Others                -0.443 0.1697    999 -2.609 
##  first_prepool Others - second_postpool Others     -0.898 0.1414    401 -6.350 
##  pool Others - second_postpool Others              -0.455 0.1251  13913 -3.638 
##  p.value
##  0.0010 
##  <.0001 
##  0.8770 
##  0.0253 
##  <.0001 
##  0.2946 
##  0.0041 
##  0.2505 
##  0.9748 
##  <.0001 
##  <.0001 
##  0.0009 
##  0.0959 
##  <.0001 
##  0.0037 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 6 estimates
summary(mod_Best)
## Linear mixed model fit by REML ['lmerMod']
## Formula: fitness ~ Step * indic_cherry + (1 | Generation)
##    Data: data_sum
## Weights: N
## 
## REML criterion at convergence: 291.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7742 -0.4669  0.1121  0.5198  2.5009 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Generation (Intercept) 0.03638  0.1907  
##  Residual               3.11306  1.7644  
## Number of obs: 247, groups:  Generation, 24
## 
## Fixed effects:
##                                        Estimate Std. Error t value
## (Intercept)                             -0.6076     0.1604  -3.788
## Steppool                                 0.9292     0.2318   4.008
## Stepsecond_postpool                      1.3086     0.1753   7.465
## indic_cherryOthers                       0.1649     0.1483   1.112
## Steppool:indic_cherryOthers             -0.4862     0.2094  -2.322
## Stepsecond_postpool:indic_cherryOthers  -0.4105     0.1605  -2.557
## 
## Correlation of Fixed Effects:
##             (Intr) Steppl Stpsc_ indc_O Stp:_O
## Steppool    -0.692                            
## Stpscnd_pst -0.915  0.633                     
## indc_chrryO -0.670  0.464  0.613              
## Stppl:ndc_O  0.475 -0.716 -0.434 -0.708       
## Stpscnd_:_O  0.619 -0.428 -0.663 -0.924  0.654

2.4 Plot with indic

#Position dodge
pd <- ggplot2::position_dodge(0.3) # move them .05 to the left and right

#Extract slope and intercept
dat_predict_allfruits <- expand.grid(Generation =as.numeric(levels(as.factor(data_sum$Generation))),
                                     Fruit_s=unique(data_sum$Fruit_s))

dat_predict_allfruits$Step <- ifelse(dat_predict_allfruits$Generation<= 7, "first_prepool",
                                     ifelse(dat_predict_allfruits$Generation>=12, "second_postpool", "pool"))
dat_predict_allfruits$indic_cherry<- ifelse(dat_predict_allfruits$Fruit_s=="Cherry", "Cherry", "Others")
dat_predict_allfruits$fitness_predicted <- predict(mod_Best, newdata = dat_predict_allfruits, 
                          re.form= NA, type = "response")

#REAL DATA
#Add G1 G6 and G7
TEMP_lineG1 <- c(rep(NA, 2), 1,rep(NA, 5))
TEMP_lineG6 <- c(rep(NA, 2), 6,rep(NA, 5))
TEMP_lineG7 <- c(rep(NA, 2), 7,rep(NA, 5))

TEMP_total <- rbind(data_sum,TEMP_lineG1,TEMP_lineG1,TEMP_lineG1,TEMP_lineG6,TEMP_lineG6,TEMP_lineG6,TEMP_lineG7,TEMP_lineG7,TEMP_lineG7)
TEMP_total$Fruit_s[TEMP_total$Generation == "6"] <- c("Strawberry", "Cranberry", "Cherry")
TEMP_total$Fruit_s[TEMP_total$Generation == "1"] <- c("Strawberry", "Cranberry", "Cherry")
TEMP_total$Fruit_s[TEMP_total$Generation == "7"] <- c("Strawberry", "Cranberry", "Cherry")
tail(TEMP_total)
##     Line    Fruit_s Generation Step  N Nb_adults sd fitness se_fitness
## 251 <NA> Strawberry          6 <NA> NA        NA NA      NA         NA
## 252 <NA>  Cranberry          6 <NA> NA        NA NA      NA         NA
## 253 <NA>     Cherry          6 <NA> NA        NA NA      NA         NA
## 254 <NA> Strawberry          7 <NA> NA        NA NA      NA         NA
## 255 <NA>  Cranberry          7 <NA> NA        NA NA      NA         NA
## 256 <NA>     Cherry          7 <NA> NA        NA NA      NA         NA
##     indic_phase1 indic_phase2 indic_phase3 indic_cherry indic_cranberry
## 251           NA            6           NA         <NA>            <NA>
## 252           NA            6           NA         <NA>            <NA>
## 253           NA            6           NA         <NA>            <NA>
## 254           NA            7           NA         <NA>            <NA>
## 255           NA            7           NA         <NA>            <NA>
## 256           NA            7           NA         <NA>            <NA>
##     indic_strawberry
## 251             <NA>
## 252             <NA>
## 253             <NA>
## 254             <NA>
## 255             <NA>
## 256             <NA>
## Add label 
TEMP_anno <- data.frame(x1 = c(3.5, 3.5, 3.5, 10, 3.5, 10), 
                        x2 = c(9, 17.5, 17.5, 17.5,  17.5, 17.5),
                        y1 = c(1, 2, 2, 1.2, 2, 1.2), 
                        y2 = c(1.25, 2.25, 2.25, 1.45,  2.25, 1.45),
                        xstar = c(6.5, 10, 10, 14,  10, 14), 
                        ystar = c(1.5, 2.5, 2.5, 1.7, 2.5, 1.7),
                        lab = c("**", "***", "***", "**",  "***", "**"),
                        Fruit_s = c("Cherry", "Cherry", 
                                    "Cranberry", "Cranberry",
                                    "Strawberry", "Strawberry"), 
                        Line = NA)
TEMP_anno
##     x1   x2  y1   y2 xstar ystar lab    Fruit_s Line
## 1  3.5  9.0 1.0 1.25   6.5   1.5  **     Cherry   NA
## 2  3.5 17.5 2.0 2.25  10.0   2.5 ***     Cherry   NA
## 3  3.5 17.5 2.0 2.25  10.0   2.5 ***  Cranberry   NA
## 4 10.0 17.5 1.2 1.45  14.0   1.7  **  Cranberry   NA
## 5  3.5 17.5 2.0 2.25  10.0   2.5 *** Strawberry   NA
## 6 10.0 17.5 1.2 1.45  14.0   1.7  ** Strawberry   NA
TEMP_title <- data.frame(xtitle = c(14, 14, 14),
                       ytitle = c(3, 3, 3),
                       title = c("Cherry", "Cranberry", "Strawberry"),
                       Fruit_s = c("Cherry", "Cranberry", "Strawberry"), 
                       Line = NA)
TEMP_title
##   xtitle ytitle      title    Fruit_s Line
## 1     14      3     Cherry     Cherry   NA
## 2     14      3  Cranberry  Cranberry   NA
## 3     14      3 Strawberry Strawberry   NA
PLOT_FITNESS_CHERRY <- ggplot(data = TEMP_total[TEMP_total$Fruit_s == "Cherry",], 
                            aes(x = factor(Generation),group = Line, y = fitness, colour =Fruit_s)) + 
  geom_errorbar(aes(ymin =fitness-1.96*se_fitness, ymax = fitness + 1.96*se_fitness),
                width=.1,position = pd, size = 0.2,color = "black") + 
  geom_line(size = 0.3,position = pd) + 
  geom_line(data = dat_predict_allfruits[dat_predict_allfruits$Fruit_s == "Cherry",], 
                     aes(x = factor(Generation), y = fitness_predicted,
                         colour = "black", group = Step), size = 0.5) + 
  geom_point(size =1, position = pd, shape =21, fill = "white") + 
  ylim(-3, 3.05) + 
  ylab("Fitness") + 
  xlab("Generation") + 
  geom_text(data = TEMP_anno[TEMP_anno$Fruit_s == "Cherry",], aes(x = xstar,  y = ystar, label = lab), size =3.3) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cherry",], aes(x = x1, xend = x1, 
           y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cherry",], aes(x = x2, xend = x2, 
           y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cherry",], aes(x = x1, xend = x2, 
           y = y2, yend = y2), size = 0.4) + 
  scale_color_manual(values = c("black", "#BC3C6D", "#FDB424", "#3FAA96")) + 
  ggtitle("Cherry") + 
  theme_LO_adaptation + theme(plot.title = element_text(color = "#BC3C6D"))
PLOT_FITNESS_CHERRY

PLOT_FITNESS_CRANB <- ggplot2::ggplot(data = TEMP_total[TEMP_total$Fruit_s == "Cranberry",], 
                            aes(x = factor(Generation),group = Line, y = fitness, colour =Fruit_s)) + 
  geom_errorbar(aes(ymin =fitness-1.96*se_fitness, ymax = fitness + 1.96*se_fitness),
                width=.1,position = pd, size = 0.2,color = "black") + 
  geom_line(size = 0.3,position = pd) + 
  geom_line(data = dat_predict_allfruits[dat_predict_allfruits$Fruit_s == "Cranberry",], aes(x = factor(Generation), y = fitness_predicted, colour = "black", group = Step),
                size = 0.5) + 
  geom_point(size =1, position = pd, shape =21, fill = "white") + 
  ylim(-3, 3.05) + 
  ylab("Fitness") + 
  xlab("Generation") + 
  geom_text(data = TEMP_anno[TEMP_anno$Fruit_s == "Cranberry",], aes(x = xstar,  y = ystar, label = lab), size =3.3) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cranberry",], aes(x = x1, xend = x1, 
           y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cranberry",], aes(x = x2, xend = x2, 
           y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cranberry",], aes(x = x1, xend = x2, 
           y = y2, yend = y2), size = 0.4) + 
  scale_color_manual(values = c("black", "#FDB424")) + 
  ggtitle("Cranberry") + 
  theme_LO_adaptation + theme(plot.title = element_text(color = "#FDB424"))
PLOT_FITNESS_CRANB

PLOT_FITNESS_STRAW <- ggplot2::ggplot(data = TEMP_total[TEMP_total$Fruit_s == "Strawberry",], 
                            aes(x = factor(Generation),group = Line, y = fitness, colour =Fruit_s)) + 
  geom_errorbar(aes(ymin =fitness-1.96*se_fitness, ymax = fitness + 1.96*se_fitness),
                width=.1,position = pd, size = 0.2,color = "black") + 
  geom_line(size = 0.3,position = pd) + 
  geom_line(data = dat_predict_allfruits[dat_predict_allfruits$Fruit_s == "Strawberry",], 
            aes(x = factor(Generation), y = fitness_predicted, colour = "black", group = Step),
            size = 0.5) + 
  geom_point(size =1, position = pd, shape =21, fill = "white") + 
  ylim(-3, 3.05) + 
  ylab("Fitness") + 
  xlab("Generation") + 
  geom_text(data = TEMP_anno[TEMP_anno$Fruit_s == "Strawberry",], 
            aes(x = xstar,  y = ystar, label = lab), size =3.3) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Strawberry",], 
               aes(x = x1, xend = x1, y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Strawberry",], 
               aes(x = x2, xend = x2, y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Strawberry",], 
               aes(x = x1, xend = x2, y = y2, yend = y2), size = 0.4) + 
  scale_color_manual(values = c("black", "#3FAA96")) + 
  ggtitle("Strawberry") + 
  theme_LO_adaptation + theme(plot.title = element_text(color = "#3FAA96"))
PLOT_FITNESS_STRAW

DYNAMIQUE_THREE_JOIN <- cowplot::ggdraw() + 
  cowplot::draw_plot(PLOT_FITNESS_CHERRY + theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank()), 
            x = 0, y = 0.66, width = 1, height = 0.33) + 
  cowplot::draw_plot(PLOT_FITNESS_CRANB + theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank()), 
            x = 0, y = 0.33, width = 1, height = 0.33) + 
  cowplot::draw_plot(PLOT_FITNESS_STRAW, 
            x = 0, y = 00, width = 1, height = 0.33) + 
  cowplot::draw_plot_label(c("A", "B", "C"),  
                  x = c(0, 0, 0), 
                  y = c(1, 0.66, 0.33), 
                  hjust = c(-0.5, -0.5, -0.5), 
                  vjust = c(1.5, 1.5, 1.5),
                  size = 12) 
 DYNAMIQUE_THREE_JOIN

cowplot::save_plot(file =here::here("figures", "FIG_TEMP_Adaptation_Indic.pdf"), DYNAMIQUE_THREE_JOIN, base_height = 17/cm(1), base_width = 11/cm(1), dpi = 1200)

3 HETEROGENEITY

3.1 Plot

pd <-  position_dodge(width = 0.5)



################### INTERMEDIATE PHENOTYPING

# Re-order levels of Line
data_sum_G7 <- data_logchange[data_logchange$Generation == "7",]
data_sum_G7$Line <- factor(data_sum_G7$Line, levels= c("CE2", "CE1", "CE4", "CE3",
                                                       "CR2", "CR3", "CR5", "CR1", "CR4",
                                                       "FR2", "FR3", "FR5", "FR1", "FR4"))


#Plot: 
symp_allop_g7 <- ggplot(data = data_sum_G7, 
                      aes(x = Line, y = logchange, group = Treatment, 
                          color = Fruit_s, shape = Treatment, fill = Fruit_s)) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size = 0.5) + 
  geom_errorbar(aes(ymin = logchange-1.96*sd_logchange, 
                    ymax =logchange + 1.96*sd_logchange, linetype = Line_type),
                  width = 0.2, size = 0.5, alpha = 0.6,position = pd) + 
  geom_point(position = pd, fill = "white", size =3) + 
  ylab("Fitness difference between intermediate\nand initial phenotyping steps")  + 
  xlab("Populations") + 
  labs(shape = "Test fruit", color = "Selection fruit") + 
  guides(color = FALSE,
           shape = guide_legend(override.aes = list(fill = c("black")))) + 
  scale_shape_manual(values = c(21, 22, 24)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
  scale_fill_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
  theme_LO_sober + theme(axis.text.x  = element_blank()) +
  coord_cartesian(expand = FALSE, ylim = c(-2, 2), xlim=c(0, 15),clip = "off") + 
  ggtitle("Intermediate phenotyping step") + 
  annotate("text", x = 2.5, y =1.5, label = 'bold(" Evolved\non cherry")',
                     size =3,colour = "#BC3C6D",parse = TRUE) + 
  annotate("text", x = 7.5, y =1.5, label = 'bold("   Evolved\non cranberry")',
                     size =3,colour = "#FDB424",parse = TRUE) + 
  annotate("text", x = 12.5, y =1.5, label = 'bold("   Evolved\non strawberry")',
                     size =3,colour = "#3FAA96",parse = TRUE) + 
  geom_segment(x = 15.5, y = -0.1, xend= 15.5, yend = -1.2, size = 0.15, 
               arrow = arrow(length = unit(0.05, "npc")),colour = "black") + 
  geom_segment(x = 15.5, y = 0.1, xend= 15.5, yend = 1.2, size = 0.15,
               arrow = arrow(length = unit(0.05, "npc")),colour = "black") + 
  annotate("text", x = 16.2, y = 0.5, label = 'bold(" Fitness\nincrease")',
             size =3,colour = "black",parse = TRUE, angle =90) + 
  annotate("text", x = 16.2, y =-0.5, label = 'bold("  Fitness\ndecrease")',
              size =3,colour = "black",parse = TRUE, angle =90) 
symp_allop_g7

################### FINAL PHENOTYPING

# Re-order levels of Line
data_sum_G29 <- data_logchange[data_logchange$Generation == "29",]
data_sum_G29$Line <- factor(data_sum_G29$Line, levels = c("CEA", "CEB", "CEC", 
                                                          "CRD", "CRA", "CRC", "CRB", "CRE",
                                                          "FRA", "FRC", "FRB"))

symp_allop_G29 <- ggplot(data = data_sum_G29, 
              aes(x = Line, y=logchange, group = Treatment, 
                  color =Fruit_s, shape = Treatment,fill =Fruit_s)) + 
      geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size = 0.5) + 
    geom_point(size =3, position = pd) + 
    geom_errorbar(aes(ymin =logchange-1.96*sd_logchange, ymax =logchange + 1.96*sd_logchange),
                width= 0.2, size = 0.5, alpha = 0.6,position = pd) + 
    ylab("Fitness difference between final\nand initial phenotyping steps")  + 
    xlab("Populations") + 
    labs(shape = "Test fruit", color = "Selection fruit") + 
    guides(fill = FALSE, alpha = FALSE) + 
    scale_shape_manual(values = c(21, 22, 24)) + 
    scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
    scale_fill_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
    theme_LO_sober + theme(axis.text.x  = element_blank()) +
    coord_cartesian(expand = FALSE, ylim = c(-2, 2), xlim=c(0, 11),clip = "off") + 
  annotate("text", x = 2, y =1.5, label = 'bold(" Evolved\non cherry")',
                     size =3,colour = "#BC3C6D",parse = TRUE) + 
  annotate("text", x = 6, y =1.5, label = 'bold("   Evolved\non cranberry")',
                     size =3,colour = "#FDB424",parse = TRUE) + 
  annotate("text", x = 10, y =1.5, label = 'bold("   Evolved\non strawberry")',
                     size =3,colour = "#3FAA96",parse = TRUE) + 
  geom_segment(x = 11.5, y = -0.1, xend= 11.5, yend = -1.2, size = 0.15, 
               arrow = arrow(length = unit(0.05, "npc")),colour = "black") + 
    geom_segment(x = 11.5, y = 0.1, xend= 11.5, yend = 1.2, size = 0.15,
               arrow = arrow(length = unit(0.05, "npc")),colour = "black") + 
    annotate("text", x = 12, y = 0.5, label = 'bold(" Fitness\nincrease")',
             size =3,colour = "black",parse = TRUE, angle =90) + 
    annotate("text", x = 12, y = -0.5, label = 'bold("  Fitness\ndecrease")',
              size =3,colour = "black", parse = TRUE, angle =90) + 
  ggtitle("Final phenotyping step") 
symp_allop_G29

### Add stroke
DIF_FITNESS_G7 <- symp_allop_g7 + 
  geom_point(aes(alpha = SA, color = interaction(SA,Fruit_s)), position = pd, size =3, stroke =1.5, fill = "white") + 
  scale_alpha_manual(values = c(0, 1)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96", "#7C2748", "#CA8702", "#328677", "#BC3C6D", "#FDB424", "#3FAA96"))  
DIF_FITNESS_G7

DIF_FITNESS_G29 <- symp_allop_G29 + geom_point(aes(alpha = SA, color =interaction(SA,Fruit_s)), position = pd, size =3, stroke =1.5) + 
      scale_alpha_manual(values = c(0, 1)) + 
      scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96", "#7C2748", "#CA8702", "#328677", "#BC3C6D", "#FDB424", "#3FAA96"))  
DIF_FITNESS_G29

plot_legend <- ggplot(data = data_logchange, aes(x = Line, y=logchange, group = Treatment, 
                  color =Fruit_s, shape = Treatment,fill =Fruit_s)) + 
    geom_point(size =3, position = pd) + 
    labs(shape = "Test fruit", color = "Selection fruit") + 
    scale_shape_manual(values = c(21, 22, 24)) + 
    scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
    theme_LO_sober + 
    guides(fill = FALSE) + 
    theme(legend.key.size = unit(0.5, "cm"),
        axis.text.x  = element_blank()) +
    coord_cartesian(expand = FALSE, ylim = c(-2, 2), xlim=c(0, 11),clip = "off") 


###########################################################
#############################################ALL PLOT
legend_trade <- lemon::g_legend(plot_legend)
    

SYMP_ALLOP_TOTAL <- cowplot::ggdraw() + 
  cowplot::draw_plot(DIF_FITNESS_G7 + theme(legend.position = 'none', 
                             plot.margin =unit(c(0.5, 3.5, 0.5, 0.5), "cm")),
            x = 0, y = 0.5, width = 1, height = 0.5) + 
  cowplot::draw_plot(DIF_FITNESS_G29 + theme(legend.position = 'none', 
                             plot.margin =unit(c(0.5, 3.5, 0.5, 0.5), "cm")), 
            x = 0, y = 0, width = 1, height = 0.5) + 
  cowplot::draw_plot_label(c("A", "B"),  
                  x = c(0, 0), y = c(1, 0.5), 
                  size = 14) + 
  cowplot::draw_plot(legend_trade, x = 0.93, y = 0.5, width = 0.001, height = 0.001) 

SYMP_ALLOP_TOTAL

cowplot::save_plot(file =here::here("figures", "FIG_Heterogeneity.pdf"),  
                   SYMP_ALLOP_TOTAL, base_height = 20/cm(1), base_width = 20/cm(1), dpi = 1200)

4 MA-AXIS REGRESSION

4.1 Analysis: a supr

# ##### Create empty dataframe with slopes and confidence interval for all pairwise
Estimates_pairwise <- data.frame(Variables = rep(rep(c("slope", "slope_CI_inf", "slope_CI_sup",
                                                          "intercept", "intercept_CI_inf", "intercept_CI_sup"),2),3),
                                    Generation = rep(c(rep(7, 6), rep(29, 6)),3),
                                    Pairwise = rep(c("Cherry_Cranberry",
                                                     "Cranberry_Strawberry", 
                                                     "Strawberry_Cherry"), each=12),
                                    Estimates = NA)


#########################
######################### G7
#########################

####### 
#######  Cherry Cranberry
####### 


mod0_G7_CheCran <- lmodel2::lmodel2(logchange_allop ~ logchange_symp,
                    range.y = "interval",range.x = "interval",
                   data = TEMP_dataG7_CheCran, nperm=1000)
mod0_G7_CheCran
## 
## Model II regression
## 
## Call: lmodel2::lmodel2(formula = logchange_allop ~ logchange_symp, data
## = TEMP_dataG7_CheCran, range.y = "interval", range.x = "interval",
## nperm = 1000)
## 
## n = 9   r = 0.52317   r-square = 0.2737069 
## Parametric P-values:   2-tailed = 0.1483644    1-tailed = 0.07418221 
## Angle between the two OLS regression lines = 30.10181 degrees
## 
## Permutation tests of OLS, MA, RMA slopes: 1-tailed, tail corresponding to sign
## A permutation test of r is equivalent to a permutation test of the OLS slope
## P-perm for SMA = NA because the SMA slope cannot be tested
## 
## Regression results
##   Method    Intercept     Slope Angle (degrees) P-perm (1-tailed)
## 1    OLS  0.015811976 0.9709274        44.15491        0.09190809
## 2     MA -0.034413319 2.8662681        70.76685        0.09190809
## 3    SMA -0.007638014 1.8558544        61.68264                NA
## 4    RMA -0.011960209 2.0189601        63.65058        0.09190809
## 
## Confidence intervals
##   Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
## 1    OLS    -0.37501980      0.40664375 -0.4426285    2.384483
## 2     MA     0.23775279      0.01609149  0.9603793   -7.404403
## 3    SMA    -0.05773723      0.01717939  0.9193256    3.746437
## 4    RMA     0.48011317      0.03818908  0.1264877  -16.550302
## 
## Eigenvalues: 0.3223501 0.05634676 
## 
## H statistic used for computing C.I. of MA: 0.2050449
#Value:
mod0_G7_CheCran$regression.results[2,]
##   Method   Intercept    Slope Angle (degrees) P-perm (1-tailed)
## 2     MA -0.03441332 2.866268        70.76685        0.09190809
mod0_G7_CheCran$confidence.intervals[2,]
##   Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
## 2     MA      0.2377528      0.01609149  0.9603793   -7.404403
#Save slope, intercept and ci
Estimates_pairwise <- filldata_maregression(model = mod0_G7_CheCran,
                                                        empty_dataset = Estimates_pairwise,
                                                        generation = 7,
                                                        pairwise = "Cherry_Cranberry")



####### 
#######  Cranberry Strawberry
####### 

mod0_G7_CranStraw <- lmodel2::lmodel2(logchange_allop ~ logchange_symp, 
                    range.y = "interval",range.x = "interval",
                   data = TEMP_dataG7_CranStraw, nperm=1000)
mod0_G7_CranStraw
## 
## Model II regression
## 
## Call: lmodel2::lmodel2(formula = logchange_allop ~ logchange_symp, data
## = TEMP_dataG7_CranStraw, range.y = "interval", range.x = "interval",
## nperm = 1000)
## 
## n = 10   r = 0.2624681   r-square = 0.0688895 
## Parametric P-values:   2-tailed = 0.4637958    1-tailed = 0.2318979 
## Angle between the two OLS regression lines = 58.18591 degrees
## 
## Permutation tests of OLS, MA, RMA slopes: 1-tailed, tail corresponding to sign
## A permutation test of r is equivalent to a permutation test of the OLS slope
## P-perm for SMA = NA because the SMA slope cannot be tested
## 
## Regression results
##   Method   Intercept     Slope Angle (degrees) P-perm (1-tailed)
## 1    OLS  0.07778912 0.4093331        22.26091         0.2337662
## 2     MA -0.62085655 3.7645194        75.12361         0.2337662
## 3    SMA -0.16171969 1.5595538        57.33164                NA
## 4    RMA -0.23805813 1.9261626        62.56315         0.2337662
## 
## Confidence intervals
##   Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
## 1    OLS     -0.2590617     0.414639991 -0.8175863    1.636253
## 2     MA      0.4815458     0.056590854  0.5111360   -1.529674
## 3    SMA     -0.5056486     0.005310832  0.7574039    3.211243
## 4    RMA             NA              NA         NA          NA
## 
## Eigenvalues: 0.09039719 0.03170788 
## 
## H statistic used for computing C.I. of MA: 0.553139
#Value: 
mod0_G7_CranStraw$confidence.intervals[2,]
##   Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
## 2     MA      0.4815458      0.05659085   0.511136   -1.529674
mod0_G7_CranStraw$regression.results[2,]
##   Method  Intercept    Slope Angle (degrees) P-perm (1-tailed)
## 2     MA -0.6208566 3.764519        75.12361         0.2337662
#Save slope, intercept and ci 
Estimates_pairwise <- filldata_maregression(model = mod0_G7_CranStraw, 
                                                        empty_dataset = Estimates_pairwise, 
                                                        generation = 7, 
                                                        pairwise = "Cranberry_Strawberry")


####### 
#######  Strawberry  Cherry
####### 
mod0_G7_StrawChe <- lmodel2::lmodel2(logchange_allop ~ logchange_symp, 
                    range.y = "interval",range.x = "interval",
                   data = TEMP_dataG7_StrawChe, nperm=1000)
mod0_G7_StrawChe
## 
## Model II regression
## 
## Call: lmodel2::lmodel2(formula = logchange_allop ~ logchange_symp, data
## = TEMP_dataG7_StrawChe, range.y = "interval", range.x = "interval",
## nperm = 1000)
## 
## n = 9   r = 0.5256197   r-square = 0.2762761 
## Parametric P-values:   2-tailed = 0.1461308    1-tailed = 0.07306542 
## Angle between the two OLS regression lines = 33.24182 degrees
## 
## Permutation tests of OLS, MA, RMA slopes: 1-tailed, tail corresponding to sign
## A permutation test of r is equivalent to a permutation test of the OLS slope
## P-perm for SMA = NA because the SMA slope cannot be tested
## 
## Regression results
##   Method Intercept     Slope Angle (degrees) P-perm (1-tailed)
## 1    OLS 0.3229824 0.3831600        20.96484        0.08091908
## 2     MA 0.3072319 0.5606481        29.27709        0.08091908
## 3    SMA 0.2922950 0.7289682        36.09086                NA
## 4    RMA 0.2854344 0.8062784        38.87848        0.08091908
## 
## Confidence intervals
##   Method 2.5%-Intercept 97.5%-Intercept  2.5%-Slope 97.5%-Slope
## 1    OLS      0.1430614       0.5029034 -0.17109307   0.9374131
## 2     MA      0.1440322       0.3707596 -0.15522634   2.3996977
## 3    SMA      0.2265349       0.3249052  0.36149308   1.4699993
## 4    RMA      0.9000034       0.3508269  0.06938925  -6.1191182
## 
## Eigenvalues: 0.1332536 0.03472534 
## 
## H statistic used for computing C.I. of MA: 0.3807415
#Value: 
mod0_G7_StrawChe$regression.results[2,]
##   Method Intercept     Slope Angle (degrees) P-perm (1-tailed)
## 2     MA 0.3072319 0.5606481        29.27709        0.08091908
mod0_G7_StrawChe$confidence.intervals[2,]
##   Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
## 2     MA      0.1440322       0.3707596 -0.1552263    2.399698
#Save slope, intercept and ci 
Estimates_pairwise <- filldata_maregression(model = mod0_G7_StrawChe, 
                                                        empty_dataset = Estimates_pairwise, 
                                                        generation = 7, 
                                                        pairwise = "Strawberry_Cherry")


#########################
######################### G29
######################### 


####### 
#######  Cherry Cranberry
####### 

mod0_G29_CheCran <- lmodel2::lmodel2(logchange_allop ~ logchange_symp, 
                    range.y = "interval",range.x = "interval",
                   data = TEMP_dataG29_CheCran, nperm=1000)
mod0_G29_CheCran
## 
## Model II regression
## 
## Call: lmodel2::lmodel2(formula = logchange_allop ~ logchange_symp, data
## = TEMP_dataG29_CheCran, range.y = "interval", range.x = "interval",
## nperm = 1000)
## 
## n = 8   r = -0.5675571   r-square = 0.322121 
## Parametric P-values:   2-tailed = 0.142274    1-tailed = 0.071137 
## Angle between the two OLS regression lines = 30.32302 degrees
## 
## Permutation tests of OLS, MA, RMA slopes: 1-tailed, tail corresponding to sign
## A permutation test of r is equivalent to a permutation test of the OLS slope
## P-perm for SMA = NA because the SMA slope cannot be tested
## 
## Regression results
##   Method Intercept      Slope Angle (degrees) P-perm (1-tailed)
## 1    OLS 0.1313351 -0.6964883       -34.85676        0.07492507
## 2     MA 0.2684057 -1.4271279       -54.98076        0.07492507
## 3    SMA 0.2308927 -1.2271686       -50.82396                NA
## 4    RMA 0.2219698 -1.1796064       -49.71071        0.07492507
## 
## Confidence intervals
##   Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
## 1    OLS    -0.20236731       0.4650375  -1.705795   0.3128181
## 2     MA     0.04073939      -1.3503555   7.201497  -0.2135780
## 3    SMA     0.10940754       0.4881066  -2.598219  -0.5796059
## 4    RMA    -0.01543750       9.6736998 -51.560993   0.0858666
## 
## Eigenvalues: 0.1689027 0.04336671 
## 
## H statistic used for computing C.I. of MA: 0.4638122
#Value: 
mod0_G29_CheCran$regression.results[2,]
##   Method Intercept     Slope Angle (degrees) P-perm (1-tailed)
## 2     MA 0.2684057 -1.427128       -54.98076        0.07492507
mod0_G29_CheCran$confidence.intervals[2,]
##   Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
## 2     MA     0.04073939       -1.350356   7.201497   -0.213578
Estimates_pairwise <- filldata_maregression(model = mod0_G29_CheCran, 
                                                        empty_dataset = Estimates_pairwise, 
                                                        generation = 29, 
                                                        pairwise = "Cherry_Cranberry")



####### 
#######  Cranberry Strawberry
####### 

mod0_G29_CranStraw <- lmodel2::lmodel2(logchange_allop ~ logchange_symp, 
                    range.y = "interval",range.x = "interval",
                   data = TEMP_dataG29_CranStraw, nperm=1000)
mod0_G29_CranStraw
## 
## Model II regression
## 
## Call: lmodel2::lmodel2(formula = logchange_allop ~ logchange_symp, data
## = TEMP_dataG29_CranStraw, range.y = "interval", range.x = "interval",
## nperm = 1000)
## 
## n = 8   r = 0.5934942   r-square = 0.3522353 
## Parametric P-values:   2-tailed = 0.1208975    1-tailed = 0.06044873 
## Angle between the two OLS regression lines = 27.62815 degrees
## 
## Permutation tests of OLS, MA, RMA slopes: 1-tailed, tail corresponding to sign
## A permutation test of r is equivalent to a permutation test of the OLS slope
## P-perm for SMA = NA because the SMA slope cannot be tested
## 
## Regression results
##   Method  Intercept     Slope Angle (degrees) P-perm (1-tailed)
## 1    OLS -0.3955687 0.7939019        38.44613        0.05794206
## 2     MA -0.5829369 1.6139086        58.21707        0.05794206
## 3    SMA -0.5198185 1.3376743        53.21945                NA
## 4    RMA -0.4703885 1.1213465        48.27390        0.05794206
## 
## Confidence intervals
##   Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
## 1    OLS     -0.8005274     0.009389961 -0.2815765    1.869380
## 2     MA      2.1350787    -0.301322822  0.3814396  -10.281344
## 3    SMA     -0.8520979    -0.360613203  0.6409208    2.791878
## 4    RMA     -1.3469425    -0.176630702 -0.1642689    4.957539
## 
## Eigenvalues: 0.2335271 0.05201112 
## 
## H statistic used for computing C.I. of MA: 0.3678656
#Value: 
mod0_G29_CranStraw$confidence.intervals[2,]
##   Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
## 2     MA       2.135079      -0.3013228  0.3814396   -10.28134
mod0_G29_CranStraw$regression.results[2,]
##   Method  Intercept    Slope Angle (degrees) P-perm (1-tailed)
## 2     MA -0.5829369 1.613909        58.21707        0.05794206
#Save slope, intercept and ci 
Estimates_pairwise <- filldata_maregression(model = mod0_G29_CranStraw, 
                                                        empty_dataset = Estimates_pairwise, 
                                                        generation = 29, 
                                                        pairwise = "Cranberry_Strawberry")


####### 
#######  Strawberry  Cherry
####### 
mod0_G29_StrawChe <- lmodel2::lmodel2(logchange_allop ~ logchange_symp, 
                    range.y = "interval",range.x = "interval",
                   data = TEMP_dataG29_StrawChe, nperm=1000)
mod0_G29_StrawChe
## 
## Model II regression
## 
## Call: lmodel2::lmodel2(formula = logchange_allop ~ logchange_symp, data
## = TEMP_dataG29_StrawChe, range.y = "interval", range.x = "interval",
## nperm = 1000)
## 
## n = 6   r = 0.6760006   r-square = 0.4569768 
## Parametric P-values:   2-tailed = 0.1404574    1-tailed = 0.07022871 
## Angle between the two OLS regression lines = 18.01702 degrees
## 
## Permutation tests of OLS, MA, RMA slopes: 1-tailed, tail corresponding to sign
## A permutation test of r is equivalent to a permutation test of the OLS slope
## P-perm for SMA = NA because the SMA slope cannot be tested
## 
## Regression results
##   Method  Intercept    Slope Angle (degrees) P-perm (1-tailed)
## 1    OLS -0.3926339 1.324565        52.94847        0.06493506
## 2     MA -0.7647082 2.537640        68.49224        0.06493506
## 3    SMA -0.5873546 1.959413        62.96220                NA
## 4    RMA -0.5806572 1.937578        62.70136        0.06393606
## 
## Confidence intervals
##   Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
## 1    OLS      -1.070322      0.28505396 -0.6798814    3.329010
## 2     MA       1.833145     -0.22100365  0.7649980   -5.932149
## 3    SMA      -1.460919     -0.23131211  0.7986067    4.807499
## 4    RMA      11.233172      0.06876396 -0.1797323  -36.579089
## 
## Eigenvalues: 0.1058706 0.01160433 
## 
## H statistic used for computing C.I. of MA: 0.2664412
#Value: 
mod0_G29_StrawChe$regression.results[2,]
##   Method  Intercept   Slope Angle (degrees) P-perm (1-tailed)
## 2     MA -0.7647082 2.53764        68.49224        0.06493506
mod0_G29_StrawChe$confidence.intervals[2,]
##   Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
## 2     MA       1.833145      -0.2210037   0.764998   -5.932149
#Save slope, intercept and ci 
Estimates_pairwise <- filldata_maregression(model = mod0_G29_StrawChe, 
                                                        empty_dataset = Estimates_pairwise, 
                                                        generation = 29, 
                                                        pairwise = "Strawberry_Cherry")


Estimates_pairwise
##           Variables Generation             Pairwise    Estimates
## 1             slope          7     Cherry_Cranberry   2.86626806
## 2      slope_CI_inf          7     Cherry_Cranberry   0.96037934
## 3      slope_CI_sup          7     Cherry_Cranberry  -7.40440288
## 4         intercept          7     Cherry_Cranberry  -0.03441332
## 5  intercept_CI_inf          7     Cherry_Cranberry   0.23775279
## 6  intercept_CI_sup          7     Cherry_Cranberry   0.01609149
## 7             slope         29     Cherry_Cranberry  -1.42712789
## 8      slope_CI_inf         29     Cherry_Cranberry   7.20149738
## 9      slope_CI_sup         29     Cherry_Cranberry  -0.21357803
## 10        intercept         29     Cherry_Cranberry   0.26840573
## 11 intercept_CI_inf         29     Cherry_Cranberry   0.04073939
## 12 intercept_CI_sup         29     Cherry_Cranberry  -1.35035553
## 13            slope          7 Cranberry_Strawberry   3.76451938
## 14     slope_CI_inf          7 Cranberry_Strawberry   0.51113600
## 15     slope_CI_sup          7 Cranberry_Strawberry  -1.52967376
## 16        intercept          7 Cranberry_Strawberry  -0.62085655
## 17 intercept_CI_inf          7 Cranberry_Strawberry   0.48154577
## 18 intercept_CI_sup          7 Cranberry_Strawberry   0.05659085
## 19            slope         29 Cranberry_Strawberry   1.61390856
## 20     slope_CI_inf         29 Cranberry_Strawberry   0.38143958
## 21     slope_CI_sup         29 Cranberry_Strawberry -10.28134362
## 22        intercept         29 Cranberry_Strawberry  -0.58293685
## 23 intercept_CI_inf         29 Cranberry_Strawberry   2.13507874
## 24 intercept_CI_sup         29 Cranberry_Strawberry  -0.30132282
## 25            slope          7    Strawberry_Cherry   0.56064812
## 26     slope_CI_inf          7    Strawberry_Cherry  -0.15522634
## 27     slope_CI_sup          7    Strawberry_Cherry   2.39969770
## 28        intercept          7    Strawberry_Cherry   0.30723191
## 29 intercept_CI_inf          7    Strawberry_Cherry   0.14403217
## 30 intercept_CI_sup          7    Strawberry_Cherry   0.37075958
## 31            slope         29    Strawberry_Cherry   2.53763976
## 32     slope_CI_inf         29    Strawberry_Cherry   0.76499799
## 33     slope_CI_sup         29    Strawberry_Cherry  -5.93214871
## 34        intercept         29    Strawberry_Cherry  -0.76470822
## 35 intercept_CI_inf         29    Strawberry_Cherry   1.83314474
## 36 intercept_CI_sup         29    Strawberry_Cherry  -0.22100365
Estimates_pairwise[Estimates_pairwise$Variables == "slope",]
##    Variables Generation             Pairwise  Estimates
## 1      slope          7     Cherry_Cranberry  2.8662681
## 7      slope         29     Cherry_Cranberry -1.4271279
## 13     slope          7 Cranberry_Strawberry  3.7645194
## 19     slope         29 Cranberry_Strawberry  1.6139086
## 25     slope          7    Strawberry_Cherry  0.5606481
## 31     slope         29    Strawberry_Cherry  2.5376398

4.2 Analysis

# ##### Create empty dataframe with slopes and confidence interval for all pairwise
Estimates_pairwise <- data.frame(Variables = rep(rep(c("slope", "slope_CI_inf", "slope_CI_sup",
                                                          "intercept"),2),3),
                                    Generation = rep(c(rep(7, 4), rep(29, 4)),3),
                                    Pairwise = rep(c("Cherry_Cranberry",
                                                     "Cranberry_Strawberry", 
                                                     "Strawberry_Cherry"), each=8),
                                    Estimates = NA)

## Nbr sample 
nb_sim = 50000


#########################
######################### G7
#########################

####### 
#######  Cherry Cranberry
####### 


# MA axis regression with bootstrap to estimate slope, intercept and 
sim <- data.frame(t(sapply(1:nb_sim, bootstrap_lmodel2, data_pairwise = TEMP_dataG7_CheCran)))

# Fill out data with quantile of simulations
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="intercept"] <- quantile(sim$intercept,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="slope"] <- quantile(sim$slope,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="slope_CI_inf"] <- quantile(sim$slope,probs = 0.025, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="slope_CI_sup"] <- quantile(sim$slope,probs = 0.975, na.rm = TRUE)


####### 
#######  Cranberry Strawberry
####### 

# MA axis regression with bootstrap to estimate slope, intercept and 
sim <- data.frame(t(sapply(1:nb_sim, bootstrap_lmodel2, data_pairwise = TEMP_dataG7_CranStraw)))

# Fill out data with quantile of simulations
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="intercept"] <- quantile(sim$intercept,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="slope"] <- quantile(sim$slope,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="slope_CI_inf"] <- quantile(sim$slope,probs = 0.025, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="slope_CI_sup"] <- quantile(sim$slope,probs = 0.975, na.rm = TRUE)


####### 
#######  Strawberry  Cherry
####### 



# MA axis regression with bootstrap to estimate slope, intercept and 
sim <- data.frame(t(sapply(1:nb_sim, bootstrap_lmodel2, data_pairwise = TEMP_dataG7_StrawChe)))

# Fill out data with quantile of simulations
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="intercept"] <- quantile(sim$intercept,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="slope"] <- quantile(sim$slope,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="slope_CI_inf"] <- quantile(sim$slope,probs = 0.025, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="slope_CI_sup"] <- quantile(sim$slope,probs = 0.975, na.rm = TRUE)


#########################
######################### G29
######################### 


####### 
#######  Cherry Cranberry
####### 

# MA axis regression with bootstrap to estimate slope, intercept and 
sim <- data.frame(t(sapply(1:nb_sim, bootstrap_lmodel2, data_pairwise = TEMP_dataG29_CheCran)))

# Fill out data with quantile of simulations
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="intercept"] <- quantile(sim$intercept,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="slope"] <- quantile(sim$slope,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="slope_CI_inf"] <- quantile(sim$slope,probs = 0.025, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="slope_CI_sup"] <- quantile(sim$slope,probs = 0.975, na.rm = TRUE)



####### 
#######  Cranberry Strawberry
####### 

# MA axis regression with bootstrap to estimate slope, intercept and 
sim <- data.frame(t(sapply(1:nb_sim, bootstrap_lmodel2, data_pairwise = TEMP_dataG29_CranStraw)))

# Fill out data with quantile of simulations
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="intercept"] <- quantile(sim$intercept,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="slope"] <- quantile(sim$slope,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="slope_CI_inf"] <- quantile(sim$slope,probs = 0.025, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="slope_CI_sup"] <- quantile(sim$slope,probs = 0.975, na.rm = TRUE)


####### 
#######  Strawberry  Cherry
####### 
rm(sim)
# MA axis regression with bootstrap to estimate slope, intercept and 
sim <- data.frame(t(sapply(1:nb_sim, bootstrap_lmodel2, data_pairwise = TEMP_dataG29_StrawChe)))

# Fill out data with quantile of simulations
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="intercept"] <- quantile(sim$intercept,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="slope"] <- quantile(sim$slope,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="slope_CI_inf"] <- quantile(sim$slope,probs = 0.025, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="slope_CI_sup"] <- quantile(sim$slope,probs = 0.975, na.rm = TRUE)



Estimates_pairwise
##       Variables Generation             Pairwise    Estimates
## 1         slope          7     Cherry_Cranberry   3.53097589
## 2  slope_CI_inf          7     Cherry_Cranberry -12.92050856
## 3  slope_CI_sup          7     Cherry_Cranberry  27.64250563
## 4     intercept          7     Cherry_Cranberry  -0.60944265
## 5         slope         29     Cherry_Cranberry  -1.63420934
## 6  slope_CI_inf         29     Cherry_Cranberry  -3.85537185
## 7  slope_CI_sup         29     Cherry_Cranberry   0.12873229
## 8     intercept         29     Cherry_Cranberry   0.24873663
## 9         slope          7 Cranberry_Strawberry   1.27496614
## 10 slope_CI_inf          7 Cranberry_Strawberry   0.41028395
## 11 slope_CI_sup          7 Cranberry_Strawberry   7.14067158
## 12    intercept          7 Cranberry_Strawberry  -0.06059237
## 13        slope         29 Cranberry_Strawberry   1.50471601
## 14 slope_CI_inf         29 Cranberry_Strawberry   0.14547472
## 15 slope_CI_sup         29 Cranberry_Strawberry   3.80483075
## 16    intercept         29 Cranberry_Strawberry  -0.49094357
## 17        slope          7    Strawberry_Cherry   1.15597786
## 18 slope_CI_inf          7    Strawberry_Cherry   0.12169649
## 19 slope_CI_sup          7    Strawberry_Cherry   3.91063251
## 20    intercept          7    Strawberry_Cherry   0.06456127
## 21        slope         29    Strawberry_Cherry   2.45170190
## 22 slope_CI_inf         29    Strawberry_Cherry   0.45700174
## 23 slope_CI_sup         29    Strawberry_Cherry   4.82393848
## 24    intercept         29    Strawberry_Cherry  -0.77826044
Estimates_pairwise[Estimates_pairwise$Variables == "slope_CI_inf"|
                     Estimates_pairwise$Variables == "slope_CI_sup",]
##       Variables Generation             Pairwise   Estimates
## 2  slope_CI_inf          7     Cherry_Cranberry -12.9205086
## 3  slope_CI_sup          7     Cherry_Cranberry  27.6425056
## 6  slope_CI_inf         29     Cherry_Cranberry  -3.8553718
## 7  slope_CI_sup         29     Cherry_Cranberry   0.1287323
## 10 slope_CI_inf          7 Cranberry_Strawberry   0.4102839
## 11 slope_CI_sup          7 Cranberry_Strawberry   7.1406716
## 14 slope_CI_inf         29 Cranberry_Strawberry   0.1454747
## 15 slope_CI_sup         29 Cranberry_Strawberry   3.8048307
## 18 slope_CI_inf          7    Strawberry_Cherry   0.1216965
## 19 slope_CI_sup          7    Strawberry_Cherry   3.9106325
## 22 slope_CI_inf         29    Strawberry_Cherry   0.4570017
## 23 slope_CI_sup         29    Strawberry_Cherry   4.8239385
########## 
########## ########## 
########## ########## 
########## ########## 
########## ########## 
########## ########## 
########## 
########## 
########## 
########## NEW FUNCTION 


temp_polygon_limits <- function(empty_dataset = Estimates_pairwise, generation = 7, pairwise = "Cherry_Cranberry") {

    # Limits for plot (arbitrary)
  ymin=-50
  ymax=50

  empty_dataset<-empty_dataset[empty_dataset$Generation==generation&
                                empty_dataset$Pairwise==pairwise,]
  
  # Create empty dataframe with slopes and confidence interval for all pairwise
  if (sign(empty_dataset$Estimates[empty_dataset$Variables=="slope_CI_inf"])!=
    sign(empty_dataset$Estimates[empty_dataset$Variables=="slope_CI_sup"])) {
      
    # Polygon limits
    intercept <- empty_dataset$Estimates[empty_dataset$Variables == "intercept"]
    inf_min<-(empty_dataset$Estimates[empty_dataset$Variables=="slope_CI_inf"]*ymin)+intercept
    inf_max<-(empty_dataset$Estimates[empty_dataset$Variables=="slope_CI_inf"]*ymax)+intercept
    sup_min <- (ymin-intercept)/empty_dataset$Estimates[empty_dataset$Variables == "slope_CI_sup"]
    sup_max <- (ymax-intercept)/empty_dataset$Estimates[empty_dataset$Variables == "slope_CI_sup"]

    # Dataframe with coordiantes
    polygon_data <- data.frame (xvalue = c(ymin, sup_min,
                                           sup_max,ymax),
                                yvalue = c(inf_min,ymin,
                                           ymax, inf_max))

  } else {
    # Polygon limits
    intercept <- empty_dataset$Estimates[empty_dataset$Variables == "intercept"]
    inf_min <- (ymin-intercept)/empty_dataset$Estimates[empty_dataset$Variables == "slope_CI_inf"]
    inf_max <- (ymax-intercept)/empty_dataset$Estimates[empty_dataset$Variables == "slope_CI_inf"]
    sup_min <- (ymin-intercept)/empty_dataset$Estimates[empty_dataset$Variables == "slope_CI_sup"]
    sup_max <- (ymax-intercept)/empty_dataset$Estimates[empty_dataset$Variables == "slope_CI_sup"]

    # Dataframe with coordiantes
    polygon_data <- data.frame (xvalue = c(inf_min,sup_min,
                                           sup_max,inf_max),
                                yvalue = c(ymin, ymin,
                                           ymax, ymax))
  }


  return(polygon_data)


}

4.3 Plot

ymin =-50
ymax =50


##################################################
##################   G7
# Calculate polygon limits and plot limits
#polygon_CheCranG7<-filldata_predict_maregression(model = mod0_G7_CheCran)
polygon_CheCranG7<-temp_polygon_limits(Estimates_pairwise, 7, "Cherry_Cranberry")

ymin_CheCranG7=min(min(TEMP_dataG7_CheCran$logchange_allop-1.96*TEMP_dataG7_CheCran$sd_allop, na.rm= TRUE),
                min(TEMP_dataG7_CheCran$logchange_symp-1.96*TEMP_dataG7_CheCran$sd_symp, na.rm= TRUE))
ymax_CheCranG7=max(max(TEMP_dataG7_CheCran$logchange_allop + 1.96*TEMP_dataG7_CheCran$sd_allop, na.rm= TRUE),
                max(TEMP_dataG7_CheCran$logchange_symp + 1.96*TEMP_dataG7_CheCran$sd_symp, na.rm= TRUE))

# Plot
CheCran_G7 <- ggplot() + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "slope"],
              colour = "black", size = 0.75) + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "slope_CI_inf"],
               colour = "grey", linetype = "dotted") + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "intercept"],
              slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "slope_CI_sup"],
              colour = "grey", linetype = "dotted") + 
  geom_errorbar(data = TEMP_dataG7_CheCran,
                aes(x = logchange_symp, ymin = logchange_allop-(1.96*sd_allop),
                    ymax = logchange_allop + (1.96*sd_allop),
                    color = Symp, linetype = Line_type),
                  width= 0.02, size = 0.3, alpha = 0.8) + 
  geom_errorbarh(data = TEMP_dataG7_CheCran,
                 aes(y = logchange_allop, xmin = logchange_symp-1.96*sd_symp, xmax = logchange_symp + 1.96*sd_symp, 
                 color = Symp, linetype = Line_type),
                 height = 0.02, size = 0.3, alpha = 0.8) + 
  geom_point(data = TEMP_dataG7_CheCran,
                 aes(x = logchange_symp, y = logchange_allop,  color = Symp,fill = Symp, shape = Allop),
                 size =3, fill = "white", stroke =1.2) + 
  xlab("Fitness difference in\nselective environment")  + 
  ylab("Fitness difference in\nalternative environment")  + 
  ggtitle("Cherry vs. Cranberry") + 
  theme_LO_sober + 
  labs(shape = "Test fruit", color = "Evolution fruit") + 
  geom_polygon(data = polygon_CheCranG7, aes(x = xvalue, y = yvalue), fill = "grey", alpha = 0.2) + 
  scale_shape_manual(values = c(21, 22)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424"))  + 
  coord_cartesian(ylim = c(ymin_CheCranG7, ymax_CheCranG7), 
                  xlim = c(ymin_CheCranG7, ymax_CheCranG7)) 
 CheCran_G7

# Calculate polygon limits and plot limits
#polygon_CranStrawG7<-filldata_predict_maregression(model = mod0_G7_CranStraw)
polygon_CranStrawG7<-temp_polygon_limits(Estimates_pairwise, 7, "Cranberry_Strawberry")

 ymin_CranStrawG7=min(min(TEMP_dataG7_CranStraw$logchange_allop-1.96*TEMP_dataG7_CranStraw$sd_allop, na.rm= TRUE),
                min(TEMP_dataG7_CranStraw$logchange_symp-1.96*TEMP_dataG7_CranStraw$sd_symp, na.rm= TRUE))
ymax_CranStrawG7=max(max(TEMP_dataG7_CranStraw$logchange_allop + 1.96*TEMP_dataG7_CranStraw$sd_allop, na.rm= TRUE),
                max(TEMP_dataG7_CranStraw$logchange_symp + 1.96*TEMP_dataG7_CranStraw$sd_symp, na.rm= TRUE))

CranStraw_G7 <-  ggplot() + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry" & 
                               Estimates_pairwise$Variables == "slope"],
              colour = "black", size = 0.75) + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry" & 
                               Estimates_pairwise$Variables == "slope_CI_inf"],
               colour = "grey", linetype = "dotted") + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry"  & 
                               Estimates_pairwise$Variables == "intercept"],
              slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7  & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry"  & 
                               Estimates_pairwise$Variables == "slope_CI_sup"],
              colour = "grey", linetype = "dotted") + 
  geom_errorbar(data = TEMP_dataG7_CranStraw,            
                aes(x =logchange_symp, ymin = logchange_allop-(1.96*sd_allop), 
                    ymax = logchange_allop + (1.96*sd_allop),
                    color = Symp, linetype = Line_type),
                  width= 0.02, size = 0.3, alpha = 0.8) + 
  geom_errorbarh(data = TEMP_dataG7_CranStraw,
                 aes(y = logchange_allop, xmin = logchange_symp-1.96*sd_symp, xmax = logchange_symp + 1.96*sd_symp, 
                 color = Symp, linetype = Line_type),
                 height = 0.02, size = 0.3, alpha = 0.8) + 
  geom_point(data = TEMP_dataG7_CranStraw,
                      aes(x =logchange_symp, y = logchange_allop,  color = Symp,fill = Symp, shape = Allop),
                   size =3, fill = "white", stroke =1.2) + 
  geom_polygon(data =polygon_CranStrawG7, aes(x = xvalue, y = yvalue), fill = "grey", alpha = 0.2) + 
  xlab("Fitness difference in\nselective environment")  + 
  ylab("Fitness difference in\nalternative environment")  + 
     ggtitle("Cranberry vs. Strawberry") + 
  coord_cartesian(ylim = c(ymin_CranStrawG7, ymax_CranStrawG7), 
                  xlim = c(ymin_CranStrawG7, ymax_CranStrawG7)) + 
   labs(shape = "Test fruit", color = "Selection fruit") + 
   scale_shape_manual(values = c(22, 24)) + 
   scale_color_manual(values = c("#FDB424", "#3FAA96"))  + 
  theme_LO_sober
 CranStraw_G7

#Calculate polygon limits and plot limits
#polygon_StrawCheG7<-filldata_predict_maregression(model = mod0_G7_StrawChe)
polygon_StrawCheG7<-temp_polygon_limits(Estimates_pairwise, 7, "Strawberry_Cherry")

 ymin_StrawCheG7=min(min(TEMP_dataG7_StrawChe$logchange_allop-1.96*TEMP_dataG7_StrawChe$sd_allop, na.rm= TRUE),
                min(TEMP_dataG7_StrawChe$logchange_symp-1.96*TEMP_dataG7_StrawChe$sd_symp, na.rm= TRUE))
ymax_StrawCheG7=max(max(TEMP_dataG7_StrawChe$logchange_allop + 1.96*TEMP_dataG7_StrawChe$sd_allop, na.rm= TRUE),
                max(TEMP_dataG7_StrawChe$logchange_symp + 1.96*TEMP_dataG7_StrawChe$sd_symp, na.rm= TRUE))

StrawChe_G7 <-  ggplot() + 
    geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "slope"],
              colour = "black", size = 0.75) + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "slope_CI_inf"],
               colour = "grey", linetype = "dotted") + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "intercept"],
              slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "slope_CI_sup"],
              colour = "grey", linetype = "dotted") + 
  geom_errorbar(data = TEMP_dataG7_StrawChe,
                aes(x =logchange_symp, ymin = logchange_allop-(1.96*sd_allop), ymax = logchange_allop + (1.96*sd_allop),
                    color = Symp, linetype = Line_type),
                  width= 0.02, size = 0.3, alpha = 0.8) + 
  geom_errorbarh(data = TEMP_dataG7_StrawChe,
                 aes(y = logchange_allop, xmin = logchange_symp-1.96*sd_symp, xmax = logchange_symp + 1.96*sd_symp, 
                 color = Symp, linetype = Line_type),
                 height = 0.02, size = 0.3, alpha = 0.8) + 
  geom_point(data = TEMP_dataG7_StrawChe,
                 aes(x =logchange_symp, y = logchange_allop,  color = Symp,fill = Symp, shape = Allop),
                 size =3, fill = "white", stroke =1.2) + 
  geom_polygon(data =polygon_StrawCheG7, aes(x = xvalue, y = yvalue), fill = "grey", alpha = 0.2) + 
    xlab("Fitness difference in\nselective environment")  + 
    ylab("Fitness difference in\nalternative environment")  + 
  ggtitle("Strawberry vs. Cherry") + 
  labs(shape = "Test fruit", color = "Selection fruit") + 
  scale_shape_manual(values = c(21, 24)) + 
  scale_color_manual(values = c("#BC3C6D", "#3FAA96"))  + 
  coord_cartesian(ylim = c(ymin_StrawCheG7, ymax_StrawCheG7), 
                  xlim = c(ymin_StrawCheG7, ymax_StrawCheG7)) + 
  theme_LO_sober
 StrawChe_G7

#####################################################################################
##################################      G29        ##################################
#####################################################################################
TEMP_dataG29_CheCran
##    Line Treatment   Fruit_s Generation   logchange sd_logchange SA      Symp
## 14  CEA Cranberry    Cherry         29 -0.21774071   0.13281291  0    Cherry
## 17  CEB Cranberry    Cherry         29 -0.59200204   0.09634940  0    Cherry
## 21  CEC Cranberry    Cherry         29 -0.20670829   0.12591289  0    Cherry
## 37  CRA    Cherry Cranberry         29  0.32537794   0.10025669  0 Cranberry
## 42  CRB    Cherry Cranberry         29  0.16929862   0.12699030  0 Cranberry
## 45  CRC    Cherry Cranberry         29  0.08474123   0.12334598  0 Cranberry
## 47  CRD    Cherry Cranberry         29  0.54736475   0.09718018  0 Cranberry
## 49  CRE    Cherry Cranberry         29 -0.10496044   0.11556266  0 Cranberry
##        Allop Line_type logchange_allop   sd_allop logchange_symp   sd_symp
## 14 Cranberry    dashed     -0.21774071 0.13281291     0.14767000 0.1574284
## 17 Cranberry    dashed     -0.59200204 0.09634940     0.25038029 0.1009864
## 21 Cranberry    dashed     -0.20670829 0.12591289     0.35854062 0.1169300
## 37    Cherry    dashed      0.32537794 0.10025669    -0.14885350 0.1193552
## 42    Cherry    dashed      0.16929862 0.12699030     0.29719689 0.1531197
## 45    Cherry    dashed      0.08474123 0.12334598     0.02703717 0.1156381
## 47    Cherry    dashed      0.54736475 0.09718018    -0.15627744 0.1156695
## 49    Cherry    dashed     -0.10496044 0.11556266     0.72513486 0.1275430
TEMP_dataG29_CranStraw
##    Line  Treatment    Fruit_s Generation   logchange sd_logchange SA       Symp
## 39  CRA Strawberry  Cranberry         29 -0.29731422    0.1098141  0  Cranberry
## 40  CRB Strawberry  Cranberry         29 -0.24539416    0.1387043  0  Cranberry
## 43  CRC Strawberry  Cranberry         29 -0.28950778    0.2012393  0  Cranberry
## 48  CRD Strawberry  Cranberry         29 -1.02555272    0.1902082  0  Cranberry
## 51  CRE Strawberry  Cranberry         29 -0.16810249    0.1398057  0  Cranberry
## 69  FRA  Cranberry Strawberry         29 -0.25158762    0.1712755  0 Strawberry
## 71  FRB  Cranberry Strawberry         29  0.04303751    0.1441119  0 Strawberry
## 74  FRC  Cranberry Strawberry         29  0.52109773    0.1254820  0 Strawberry
##         Allop Line_type logchange_allop  sd_allop logchange_symp   sd_symp
## 39 Strawberry    dashed     -0.29731422 0.1098141    -0.14885350 0.1193552
## 40 Strawberry    dashed     -0.24539416 0.1387043     0.29719689 0.1531197
## 43 Strawberry    dashed     -0.28950778 0.2012393     0.02703717 0.1156381
## 48 Strawberry    dashed     -1.02555272 0.1902082    -0.15627744 0.1156695
## 51 Strawberry    dashed     -0.16810249 0.1398057     0.72513486 0.1275430
## 69  Cranberry    dashed     -0.25158762 0.1712755     0.15709032 0.1411560
## 71  Cranberry    dashed      0.04303751 0.1441119     0.56023100 0.1284283
## 74  Cranberry    dashed      0.52109773 0.1254820     0.36640738 0.1201170
TEMP_dataG29_StrawChe
##    Line  Treatment    Fruit_s Generation   logchange sd_logchange SA       Symp
## 13  CEA Strawberry     Cherry         29  0.01971359    0.1475176  0     Cherry
## 16  CEB Strawberry     Cherry         29 -0.15850680    0.0983504  0     Cherry
## 19  CEC Strawberry     Cherry         29  0.21131464    0.1237821  0     Cherry
## 68  FRA     Cherry Strawberry         29 -0.50599095    0.1378652  0 Strawberry
## 72  FRB     Cherry Strawberry         29  0.19377147    0.1094819  0 Strawberry
## 75  FRC     Cherry Strawberry         29  0.32151693    0.1140574  0 Strawberry
##         Allop Line_type logchange_allop  sd_allop logchange_symp   sd_symp
## 13 Strawberry    dashed      0.01971359 0.1475176      0.1476700 0.1574284
## 16 Strawberry    dashed     -0.15850680 0.0983504      0.2503803 0.1009864
## 19 Strawberry    dashed      0.21131464 0.1237821      0.3585406 0.1169300
## 68     Cherry    dashed     -0.50599095 0.1378652      0.1570903 0.1411560
## 72     Cherry    dashed      0.19377147 0.1094819      0.5602310 0.1284283
## 75     Cherry    dashed      0.32151693 0.1140574      0.3664074 0.1201170
#Calculate polygon limits and plot limits
polygon_CheCranG29<-filldata_predict_maregression(model = mod0_G29_CheCran)
polygon_CheCranG29<-temp_polygon_limits(Estimates_pairwise, 29, "Cherry_Cranberry")

ymin_CheCranG29=min(min(TEMP_dataG29_CheCran$logchange_allop-1.96*TEMP_dataG29_CheCran$sd_allop, na.rm= TRUE),
                min(TEMP_dataG29_CheCran$logchange_symp-1.96*TEMP_dataG29_CheCran$sd_symp, na.rm= TRUE))
ymax_CheCranG29=max(max(TEMP_dataG29_CheCran$logchange_allop + 1.96*TEMP_dataG29_CheCran$sd_allop, na.rm= TRUE),
                max(TEMP_dataG29_CheCran$logchange_symp + 1.96*TEMP_dataG29_CheCran$sd_symp, na.rm= TRUE))


CheCran_G29 <- ggplot() + 
    geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "slope"],
              colour = "black", size = 0.75) + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "slope_CI_inf"],
               colour = "grey", linetype = "dotted") + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "intercept"],
              slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "slope_CI_sup"],
              colour = "grey", linetype = "dotted") + 
  geom_errorbar(data = TEMP_dataG29_CheCran,
                aes(x =logchange_symp, ymin = logchange_allop-(1.96*sd_allop), 
                    ymax = logchange_allop + (1.96*sd_allop),
                    color = Symp),
                  width= 0.02, size = 0.3, alpha = 0.8) + 
  geom_errorbarh(data = TEMP_dataG29_CheCran,
                 aes(y = logchange_allop, xmin = logchange_symp-1.96*sd_symp, xmax = logchange_symp + 1.96*sd_symp, 
                 color = Symp),
                 height = 0.02, size = 0.3, alpha = 0.8) + 
  geom_point(data = TEMP_dataG29_CheCran,
                 aes(x =logchange_symp, y = logchange_allop,  color = Symp,fill = Symp, shape = Allop),
                 size =3, fill = "white", stroke =1.2) + 
  geom_polygon(data =polygon_CheCranG29, aes(x = xvalue, y = yvalue), fill = "grey", alpha = 0.2) + 
    xlab("Fitness difference in\nselective environment")  + 
    ylab("Fitness difference in\nalternative environment")  + 
  ggtitle("Cherry vs. Cranberry") + 
  labs(shape = "Test fruit", color = "Evolution fruit") + 
  scale_shape_manual(values = c(16, 15)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424"))  + 
  coord_cartesian(ylim = c(ymin_CheCranG29, ymax_CheCranG29), 
                  xlim = c(ymin_CheCranG29, ymax_CheCranG29)) + 
  theme_LO_sober
 CheCran_G29

#Calculate polygon limits and plot limits
polygon_CranStrawG29<-filldata_predict_maregression(model = mod0_G29_CranStraw)
polygon_CranStrawG29<-temp_polygon_limits(Estimates_pairwise, 29, "Cranberry_Strawberry")

ymin_CranStrawG29=min(min(TEMP_dataG29_CranStraw$logchange_allop-1.96*TEMP_dataG29_CranStraw$sd_allop, na.rm= TRUE),
                min(TEMP_dataG29_CranStraw$logchange_symp-1.96*TEMP_dataG29_CranStraw$sd_symp, na.rm= TRUE))
ymax_CranStrawG29=max(max(TEMP_dataG29_CranStraw$logchange_allop + 1.96*TEMP_dataG29_CranStraw$sd_allop, na.rm= TRUE),
                max(TEMP_dataG29_CranStraw$logchange_symp + 1.96*TEMP_dataG29_CranStraw$sd_symp, na.rm= TRUE))

#Plot
 CranStraw_G29 <- ggplot() + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry" & 
                               Estimates_pairwise$Variables == "slope"],
              colour = "black", size = 0.75) + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry" & 
                               Estimates_pairwise$Variables == "slope_CI_inf"],
               colour = "grey", linetype = "dotted") + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry" & 
                               Estimates_pairwise$Variables == "intercept"],
              slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry" & 
                               Estimates_pairwise$Variables == "slope_CI_sup"],
                            colour = "grey", linetype = "dotted") + 
   geom_errorbar(data = TEMP_dataG29_CranStraw,
                aes(x =logchange_symp, ymin = logchange_allop-(1.96*sd_allop), 
                    ymax = logchange_allop + (1.96*sd_allop),
                    color = Symp),
                  width= 0.02, size = 0.3, alpha = 0.8) + 
   geom_errorbarh(data = TEMP_dataG29_CranStraw,
                 aes(y = logchange_allop, xmin = logchange_symp-1.96*sd_symp, xmax = logchange_symp + 1.96*sd_symp, 
                 color = Symp),
                 height = 0.02, size = 0.3, alpha = 0.8) + 
   geom_point(data = TEMP_dataG29_CranStraw,
                 aes(x =logchange_symp, y = logchange_allop,  color = Symp,fill = Symp, shape = Allop),
                 size =3, fill = "white", stroke =1.2) + 
   geom_polygon(data =polygon_CranStrawG29, aes(x = xvalue, y = yvalue), fill = "grey", alpha = 0.2) + 
    xlab("Fitness difference in\nselective environment")  + 
    ylab("Fitness difference in\nalternative environment")  + 
   ggtitle("Cranberry vs. Strawberry") + 
   coord_cartesian(ylim = c(ymin_CranStrawG29, ymax_CranStrawG29), 
                  xlim = c(ymin_CranStrawG29, ymax_CranStrawG29)) + 
   labs(shape = "Test fruit", color = "Selection fruit") + 
   scale_shape_manual(values = c(15, 17)) + 
   scale_color_manual(values = c("#FDB424", "#3FAA96"))  + 
   theme_LO_sober
 CranStraw_G29

#Calculate polygon limits and plot limits
polygon_StrawCheG29<-filldata_predict_maregression(model = mod0_G29_StrawChe)
polygon_StrawCheG29<-temp_polygon_limits(Estimates_pairwise, 29, "Strawberry_Cherry")

ymin_StrawCheG29=min(min(TEMP_dataG29_StrawChe$logchange_allop-1.96*TEMP_dataG29_StrawChe$sd_allop, na.rm= TRUE),
                min(TEMP_dataG29_StrawChe$logchange_symp-1.96*TEMP_dataG29_StrawChe$sd_symp, na.rm= TRUE))
ymax_StrawCheG29=max(max(TEMP_dataG29_StrawChe$logchange_allop + 1.96*TEMP_dataG29_StrawChe$sd_allop, na.rm= TRUE),
                max(TEMP_dataG29_StrawChe$logchange_symp + 1.96*TEMP_dataG29_StrawChe$sd_symp, na.rm= TRUE))


StrawChe_G29 <- ggplot() + 
 geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "slope"],
              colour = "black", size = 0.75) + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "slope_CI_inf"],
               colour = "grey", linetype = "dotted") + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "intercept"],
              slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "slope_CI_sup"],
                            colour = "grey", linetype = "dotted") + 
  geom_errorbar(data = TEMP_dataG29_StrawChe,
                aes(x =logchange_symp, ymin = logchange_allop-(1.96*sd_allop), 
                    ymax = logchange_allop + (1.96*sd_allop),
                    color = Symp),
                  width= 0.02, size = 0.2, alpha =1) + 
   geom_errorbarh(data = TEMP_dataG29_StrawChe,
                 aes(y = logchange_allop, xmin = logchange_symp-1.96*sd_symp, xmax = logchange_symp + 1.96*sd_symp, 
                 color = Symp),
                 height = 0.02, size = 0.2, alpha =1) + 
  geom_point(data = TEMP_dataG29_StrawChe,
                 aes(x =logchange_symp, y = logchange_allop,  color = Symp,fill = Symp, shape = Allop),
                 size =3, fill = "white", stroke =1.2) + 
  geom_polygon(data =polygon_StrawCheG29, aes(x = xvalue, y = yvalue), fill = "grey", alpha = 0.2) + 
  coord_cartesian(ylim = c(ymin_StrawCheG29, ymax_StrawCheG29), 
                  xlim = c(ymin_StrawCheG29, ymax_StrawCheG29)) + 
    xlab("Fitness difference in\nselective environment")  + 
    ylab("Fitness difference in\nalternative environment")  + 
  ggtitle("Strawberry vs. Cherry") + 
  labs(shape = "Test fruit", color = "Selection fruit") + 
  scale_shape_manual(values = c(16, 17)) + 
  scale_color_manual(values = c("#BC3C6D", "#3FAA96"))  + 
  theme_LO_sober
 StrawChe_G29

legend_tradevide <-  ggplot(data = data_logchange[data_logchange$Generation == "29",],
                          aes(x =logchange, y = logchange,  color = Symp,fill = Symp, shape = Allop)) + 
  geom_point(size =2.5, fill = "white") + 
  labs(shape = "Test fruit", color = "Selection fruit") + 
  scale_shape_manual(values = c(16, 15, 17)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96"))  + 
  theme_LO_sober

legend_trade <- lemon::g_legend(legend_tradevide)

 
 

Slopeestimates_logchange_pairwise_errorbar <- cowplot::ggdraw() + 
  cowplot::draw_plot(CheCran_G7 + theme(legend.position = "none"), 
            x = 0.00, y = 0.5, width = 0.22, height = 0.45) + 
  cowplot::draw_plot(CranStraw_G7 + theme(legend.position = "none"), 
            x = 0.26, y = 0.5, width = 0.22, height = 0.45) + 
  cowplot::draw_plot(StrawChe_G7 + theme(legend.position = "none"), 
            x = 0.52, y = 0.5, width = 0.22, height = 0.45) + 
  cowplot::draw_plot(legend_trade, x = 0.750, y = 0.5, width = 0.1, height = 0.1) + 
  cowplot::draw_plot(CheCran_G29 + theme(legend.position = "none"), 
            x = 0.00, y = 0, width = 0.22, height = 0.45) + 
  cowplot::draw_plot(CranStraw_G29 + theme(legend.position = "none"), 
            x = 0.26, y = 0, width = 0.22, height = 0.45) + 
  cowplot::draw_plot(StrawChe_G29 + theme(legend.position = "none"), 
            x = 0.52, y = 0, width = 0.22, height = 0.45) + 
  cowplot::draw_plot_label(c("Intermediate phenotyping step", "A", "B", "C", " ",
                    "Final phenotyping step", "D", "E", "F", " "),  
                  x = c(0.26, 0, 0.26, 0.52, 0.82, 0.2, 0, 0.26, 0.52, 0.82), 
                  y = c(1, 0.95, 0.95, 0.95, 0.95, 0.5, 0.45, 0.45, 0.45, 0.45), 
                  hjust = c(-0.25, -0.25, -0.25, -0.25, -0.25, -0.75, -0.75, -0.75, -0.75, -0.75), 
                  vjust = c(1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5),
                  size = 14) 
 Slopeestimates_logchange_pairwise_errorbar

#  
# #  
# 
cowplot::save_plot(file =here::here("figures", "FIG_MAregression.pdf"),
                   Slopeestimates_logchange_pairwise_errorbar,
                   base_height = 15/cm(1), base_width = 30/cm(1), dpi = 610)

5 SA AND LOCAL ADAPTATION

5.1 Analysis

##### 
#####  G7
##### 
# Pblm when Nb_adults= 0, it is impossible to calculate fitness log(0/20)
data_G7[data_G7$Nb_adults == 0,]
##      Treatment Line    Fruit_s Nb_eggs Nb_adults SA Emergence_rate
## 85  Strawberry  FR4 Strawberry      28         0  1              0
## 158  Cranberry  CE1     Cherry     365         0  0              0
## 159     Cherry  FR1 Strawberry     284         0  0              0
## 177     Cherry  CE1     Cherry     249         0  1              0
# correspond to 2 allopatric and 2 sympatric measurements (2 in CE1, 1 in FR4 et 1 in FR1)

# Remove Nb_adults = 0 because it is negligeable
data_G7_without0 <- data_G7[data_G7$Nb_adults!= 0,]
data_G7_without0$fitness <- log(data_G7_without0$Nb_adults/20)

lm_val_G7 = lm(fitness ~ Treatment + Line + SA + Treatment:Fruit_s, data = data_G7_without0)

Fratio = anova(lm_val_G7)[3, 3]/anova(lm_val_G7)[4, 3]
pvalue = 1 - pf(Fratio, anova(lm_val_G7)[3, 1], anova(lm_val_G7)[4, 1]) 

pvalue
## [1] 0.9570523
Fratio
## [1] 0.003418531
anova(lm_val_G7)[3, 1]
## [1] 1
anova(lm_val_G7)[4, 1]
## [1] 3
##### 
#####  G29
##### 

# Pblm when Nb_adults= 0, it is impossible to calculate fitness log(0/20)
data_G29[data_G29$Nb_adults == 0,]
##       Treatment Line    Fruit_s Nb_eggs Nb_adults SA Emergence_rate
## 5406     Cherry  CEA     Cherry      18         0  1              0
## 5494 Strawberry  CRA  Cranberry      56         0  0              0
## 5501 Strawberry  CRA  Cranberry      82         0  0              0
## 5535 Strawberry  CEB     Cherry      74         0  0              0
## 5555     Cherry  CRE  Cranberry      19         0  0              0
## 5563 Strawberry  CRE  Cranberry     113         0  0              0
## 5564 Strawberry  CRE  Cranberry     125         0  0              0
## 5629 Strawberry  CRB  Cranberry     106         0  0              0
## 5741     Cherry  CEA     Cherry      23         0  1              0
## 5769  Cranberry  FRB Strawberry      23         0  0              0
## 5831  Cranberry  CRA  Cranberry      40         0  1              0
## 5922 Strawberry  FRA Strawberry      82         0  1              0
## 5997 Strawberry  CRD  Cranberry      25         0  0              0
## 6059     Cherry  CEA     Cherry     187         0  1              0
## 6062 Strawberry  CEA     Cherry     195         0  0              0
## 6090 Strawberry  FRB Strawberry     143         0  1              0
## 6111     Cherry  FRB Strawberry     114         0  0              0
## 6240  Cranberry  FRA Strawberry     125         0  0              0
## 6293     Cherry  CEC     Cherry     188         0  1              0
## 6301     Cherry  CEC     Cherry     163         0  1              0
## 6309  Cranberry  CEC     Cherry     199         0  0              0
## 6324 Strawberry  CRD  Cranberry      82         0  0              0
## 6325 Strawberry  CRD  Cranberry      72         0  0              0
## 6327 Strawberry  CRD  Cranberry      58         0  0              0
## 6330 Strawberry  CRD  Cranberry      59         0  0              0
length(data_G29[data_G29$Nb_adults == 0 & data_G29$SA == 0,]$Nb_adults)
## [1] 17
length(data_G29[data_G29$Nb_adults == 0 & data_G29$SA == 1,]$Nb_adults)
## [1] 8
# correspond to 25 measures 17 allopatric and 8 sympatric measurements (2 in CE1, 1 in FR4 et 1 in FR1)

# Remove Nb_adults = 0 because it is negligeable
data_G29_without0 <- data_G29[data_G29$Nb_adults!= 0,]
data_G29_without0$fitness <- log(data_G29_without0$Nb_adults/20)

lm_val_G29 = lm(fitness ~ Treatment + Line + SA + Treatment:Fruit_s, data = data_G29_without0)

Fratio = anova(lm_val_G29)[3, 3]/anova(lm_val_G29)[4, 3]
pvalue = 1 - pf(Fratio, anova(lm_val_G29)[3, 1], anova(lm_val_G29)[4, 1]) 

pvalue
## [1] 0.1733915
Fratio
## [1] 3.1628
anova(lm_val_G29)[3, 1]
## [1] 1
anova(lm_val_G29)[4, 1]
## [1] 3